From 6382eea05bf1d8abaa18c0397c80a4982866e575 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Thu, 5 Dec 2024 19:43:05 -0500 Subject: [PATCH 01/12] Add documentation to the functions. Signed-off-by: Michael Jackson --- .../OrientationAnalysis/CMakeLists.txt | 4 +- .../Filters/Algorithms/ComputeShapes.cpp | 39 ++-- .../Filters/Algorithms/ComputeShapes.hpp | 12 +- .../Algorithms/ComputeTriangleGeomShapes.cpp | 172 +++++++++++------- 4 files changed, 132 insertions(+), 95 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/CMakeLists.txt b/src/Plugins/OrientationAnalysis/CMakeLists.txt index 4d77daf939..adf68735bf 100644 --- a/src/Plugins/OrientationAnalysis/CMakeLists.txt +++ b/src/Plugins/OrientationAnalysis/CMakeLists.txt @@ -52,7 +52,7 @@ set(FilterList ComputeSchmidsFilter ComputeShapesFilter ComputeSlipTransmissionMetricsFilter - # ComputeTriangleGeomShapesFilter + ComputeTriangleGeomShapesFilter ConvertHexGridToSquareGridFilter ConvertOrientationsFilter ConvertQuaternionFilter @@ -186,7 +186,7 @@ set(filter_algorithms ComputeSchmids ComputeShapes ComputeSlipTransmissionMetrics - # ComputeTriangleGeomShapes + ComputeTriangleGeomShapes ConvertHexGridToSquareGrid ConvertQuaternion CreateEnsembleInfo diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.cpp index 834dd39a9d..773e3e1fca 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.cpp @@ -139,27 +139,28 @@ Result<> ComputeShapes::operator()() if(imageGeom.getNumXCells() > 1 && imageGeom.getNumYCells() > 1 && imageGeom.getNumZCells() > 1) { - find_moments(); - find_axes(); - find_axiseulers(); + findMoments(); + findAxes(); + findAxisEulers(); } if(imageGeom.getNumXCells() == 1 || imageGeom.getNumYCells() == 1 || imageGeom.getNumZCells() == 1) { - find_moments2D(); - find_axes2D(); - find_axiseulers2D(); + findMoments2D(); + findAxes2D(); + findAxisEulers2D(); } return {}; } // ----------------------------------------------------------------------------- -void ComputeShapes::find_moments() +void ComputeShapes::findMoments() { const auto& imageGeom = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryPath); const auto& featureIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureIdsArrayPath); const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); + // Calculated Arrays auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); auto& omega3s = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); @@ -208,9 +209,11 @@ void ComputeShapes::find_moments() y2 = y - (modYRes / 4.0f); z1 = z + (modZRes / 4.0f); z2 = z - (modZRes / 4.0f); + xdist1 = (x1 - (centroids[gnum * 3 + 0] * static_cast(m_ScaleFactor))); ydist1 = (y1 - (centroids[gnum * 3 + 1] * static_cast(m_ScaleFactor))); zdist1 = (z1 - (centroids[gnum * 3 + 2] * static_cast(m_ScaleFactor))); + xdist2 = (x1 - (centroids[gnum * 3 + 0] * static_cast(m_ScaleFactor))); ydist2 = (y1 - (centroids[gnum * 3 + 1] * static_cast(m_ScaleFactor))); zdist2 = (z2 - (centroids[gnum * 3 + 2] * static_cast(m_ScaleFactor))); @@ -296,6 +299,7 @@ void ComputeShapes::find_moments() m_FeatureEigenVals[featureId * 3 + 1] = eigenValues[idxs[1]].real(); m_FeatureEigenVals[featureId * 3 + 2] = eigenValues[idxs[2]].real(); + // These values will be used to compute the axis eulers // EigenVector associated with the largest EigenValue goes in the 3rd column auto col = eigenVectors.col(idxs[0]); m_EFVec[featureId * 9 + 2] = col(0).real(); @@ -321,6 +325,7 @@ void ComputeShapes::find_moments() u110 = static_cast(-m_FeatureMoments[featureId * 6 + 3]); u011 = static_cast(-m_FeatureMoments[featureId * 6 + 4]); u101 = static_cast(-m_FeatureMoments[featureId * 6 + 5]); + o3 = static_cast((u200 * u020 * u002) + (2.0f * u110 * u101 * u011) - (u200 * u011 * u011) - (u020 * u101 * u101) - (u002 * u110 * u110)); vol5 = pow(vol5, 5.0); omega3 = vol5 / o3; @@ -338,9 +343,7 @@ void ComputeShapes::find_moments() } // ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- -void ComputeShapes::find_moments2D() +void ComputeShapes::findMoments2D() { const auto& featureIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureIdsArrayPath); @@ -430,9 +433,7 @@ void ComputeShapes::find_moments2D() } // ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- -void ComputeShapes::find_axes() +void ComputeShapes::findAxes() { auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); auto& aspectRatios = m_DataStructure.getDataRefAs(m_InputValues->AspectRatiosArrayPath); @@ -476,9 +477,7 @@ void ComputeShapes::find_axes() } // ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- -void ComputeShapes::find_axes2D() +void ComputeShapes::findAxes2D() { auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); @@ -551,9 +550,7 @@ void ComputeShapes::find_axes2D() } // ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- -void ComputeShapes::find_axiseulers() +void ComputeShapes::findAxisEulers() { const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); auto& axisEulerAngles = m_DataStructure.getDataRefAs(m_InputValues->AxisEulerAnglesArrayPath); @@ -588,9 +585,7 @@ void ComputeShapes::find_axiseulers() } // ----------------------------------------------------------------------------- -// -// ----------------------------------------------------------------------------- -void ComputeShapes::find_axiseulers2D() +void ComputeShapes::findAxisEulers2D() { const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); auto& axisEulerAngles = m_DataStructure.getDataRefAs(m_InputValues->AxisEulerAnglesArrayPath); diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.hpp index d425b4f837..22c4a783a0 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapes.hpp @@ -56,32 +56,32 @@ class ORIENTATIONANALYSIS_EXPORT ComputeShapes /** * @brief find_moments Determines the second order moments for each Feature */ - void find_moments(); + void findMoments(); /** * @brief find_moments2D Determines the second order moments for each Feature (2D version) */ - void find_moments2D(); + void findMoments2D(); /** * @brief find_axes Determine principal axis lengths for each Feature */ - void find_axes(); + void findAxes(); /** * @brief find_axes2D Determine principal axis lengths for each Feature (2D version) */ - void find_axes2D(); + void findAxes2D(); /** * @brief find_axiseulers Determine principal axis directions for each Feature */ - void find_axiseulers(); + void findAxisEulers(); /** * @brief find_axiseulers2D Determine principal axis directions for each Feature (2D version) */ - void find_axiseulers2D(); + void findAxisEulers2D(); }; } // namespace nx::core diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 9211f382b5..a155ec6524 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -16,50 +16,46 @@ constexpr float64 k_ScaleFactor = 1.0; constexpr usize k_00 = 0; constexpr usize k_01 = 1; constexpr usize k_02 = 2; + constexpr usize k_10 = 3; constexpr usize k_11 = 4; constexpr usize k_12 = 5; + constexpr usize k_20 = 6; constexpr usize k_21 = 7; constexpr usize k_22 = 8; // ----------------------------------------------------------------------------- +/** + * @brief This function will calculate the volume and centroid for 8 sub tetrahedrons that are generated + * from the triangle given by the "verts" and a centroid point give by "centroid" + * @tparam T + * @param verts The 3 Vertices of the triangle + * @param centroid Another point in space that represents the centroid of the feature that the triangle belongs to + * @param tetInfo The array of TetInfo data which is the volume and centroid of each sub-tetrahedron that is formed. + */ template -void FindTetrahedronInfo(const std::array& vertIds, const DataArray& vertPtr, const std::array& centroid, std::array& tetInfo) +void FindTetrahedronInfo(const std::array& verts, const std::array& centroid, std::array& tetInfo) { - std::array coords = {vertPtr[3 * vertIds[0] + 0], - vertPtr[3 * vertIds[0] + 1], - vertPtr[3 * vertIds[0] + 2], - vertPtr[3 * vertIds[1] + 0], - vertPtr[3 * vertIds[1] + 1], - vertPtr[3 * vertIds[1] + 2], - vertPtr[3 * vertIds[2] + 0], - vertPtr[3 * vertIds[2] + 1], - vertPtr[3 * vertIds[2] + 2], - centroid[0], - centroid[1], - centroid[2], - 0.5 * (vertPtr[3 * vertIds[0] + 0] + centroid[0]), - 0.5 * (vertPtr[3 * vertIds[0] + 1] + centroid[1]), - 0.5 * (vertPtr[3 * vertIds[0] + 2] + centroid[2]), - 0.5 * (vertPtr[3 * vertIds[1] + 0] + centroid[0]), - 0.5 * (vertPtr[3 * vertIds[1] + 1] + centroid[1]), - 0.5 * (vertPtr[3 * vertIds[1] + 2] + centroid[2]), - 0.5 * (vertPtr[3 * vertIds[2] + 0] + centroid[0]), - 0.5 * (vertPtr[3 * vertIds[2] + 1] + centroid[1]), - 0.5 * (vertPtr[3 * vertIds[2] + 2] + centroid[2]), - 0.5 * (vertPtr[3 * vertIds[0] + 0] + vertPtr[3 * vertIds[1] + 0]), - 0.5 * (vertPtr[3 * vertIds[0] + 1] + vertPtr[3 * vertIds[1] + 1]), - 0.5 * (vertPtr[3 * vertIds[0] + 2] + vertPtr[3 * vertIds[1] + 2]), - 0.5 * (vertPtr[3 * vertIds[1] + 0] + vertPtr[3 * vertIds[2] + 0]), - 0.5 * (vertPtr[3 * vertIds[1] + 1] + vertPtr[3 * vertIds[2] + 1]), - 0.5 * (vertPtr[3 * vertIds[1] + 2] + vertPtr[3 * vertIds[2] + 2]), - 0.5 * (vertPtr[3 * vertIds[2] + 0] + vertPtr[3 * vertIds[0] + 0]), - 0.5 * (vertPtr[3 * vertIds[2] + 1] + vertPtr[3 * vertIds[0] + 1]), - 0.5 * (vertPtr[3 * vertIds[2] + 2] + vertPtr[3 * vertIds[0] + 2])}; // clang-format off + std::array coords = {verts[0][0], verts[0][1], verts[0][2], // V0 + verts[1][0], verts[1][1], verts[1][2], // V1 + verts[2][0], verts[2][1], verts[2][2], // V2 + + centroid[0], centroid[1], centroid[2], // Centroid Vertex + + 0.5 * (verts[0][0] + centroid[0]), 0.5 * (verts[0][1] + centroid[1]), 0.5 * (verts[0][2] + centroid[2]), // Average of V0 and Centroid + 0.5 * (verts[1][0] + centroid[0]), 0.5 * (verts[1][1] + centroid[1]), 0.5 * (verts[1][2] + centroid[2]), // Average of V1 and Centroid + 0.5 * (verts[2][0] + centroid[0]), 0.5 * (verts[2][1] + centroid[1]), 0.5 * (verts[2][2] + centroid[2]), // Average of V2 and Centroid + + 0.5 * (verts[0][0] + verts[1][0]), 0.5 * (verts[0][1] + verts[1][1]), 0.5 * (verts[0][2] + verts[1][2]), // Average of V0 & V1 Distance + 0.5 * (verts[1][0] + verts[2][0]), 0.5 * (verts[1][1] + verts[2][1]), 0.5 * (verts[1][2] + verts[2][2]), // Average of V1 & V2 Distance + 0.5 * (verts[2][0] + verts[0][0]), 0.5 * (verts[2][1] + verts[0][1]), 0.5 * (verts[2][2] + verts[0][2]) // Average of V2 & V0 Distance + }; + + // Create 8 Tetrahedrons std::array tets = { - 4, 5, 6, 3, + 4, 5, 6, 3, 0, 7, 9, 4, 1, 8, 7, 5, 2, 9, 8, 6, @@ -70,48 +66,77 @@ void FindTetrahedronInfo(const std::array& vertIds, const DataArray }; // clang-format on + using Point3DType = Point3D; + // Loop over each Tetrahedron for(usize iter = 0; iter < 8; iter++) { - T ax = coords[3 * tets[4 * iter + 0] + 0]; - T ay = coords[3 * tets[4 * iter + 0] + 1]; - T az = coords[3 * tets[4 * iter + 0] + 2]; - T bx = coords[3 * tets[4 * iter + 1] + 0]; - T by = coords[3 * tets[4 * iter + 1] + 1]; - T bz = coords[3 * tets[4 * iter + 1] + 2]; - T cx = coords[3 * tets[4 * iter + 2] + 0]; - T cy = coords[3 * tets[4 * iter + 2] + 1]; - T cz = coords[3 * tets[4 * iter + 2] + 2]; - T dx = coords[3 * tets[4 * iter + 3] + 0]; - T dy = coords[3 * tets[4 * iter + 3] + 1]; - T dz = coords[3 * tets[4 * iter + 3] + 2]; - - std::array volumeMatrix = {bx - ax, cx - ax, dx - ax, by - ay, cy - ay, dy - ay, bz - az, cz - az, dz - az}; - + usize i0 = 4 * iter + 0; + usize i1 = 4 * iter + 1; + usize i2 = 4 * iter + 2; + usize i3 = 4 * iter + 3; + // Create the 4 Vertices of the sub-tetrahedron + Point3DType a(coords[3 * tets[i0] + 0], coords[3 * tets[i0] + 1], coords[3 * tets[i0] + 2]); + Point3DType b(coords[3 * tets[i1] + 0], coords[3 * tets[i1] + 1], coords[3 * tets[i1] + 2]); + Point3DType c(coords[3 * tets[i2] + 0], coords[3 * tets[i2] + 1], coords[3 * tets[i2] + 2]); + Point3DType d(coords[3 * tets[i3] + 0], coords[3 * tets[i3] + 1], coords[3 * tets[i3] + 2]); + + // The volume of the tetrahedron can be calculated through V= 1/6 |M| (where M is the volume matrix) + // v1 = b-a, v2 = c-a, v3 = d-a + // | v1x v2x v3x + // M = | v1y v2y v3y + // | v1z v2z v3z + + // clang-format off + std::array volumeMatrix = {b[0] - a[0], c[0] - a[0], d[0] - a[0], + b[1] - a[1], c[1] - a[1], d[1] - a[1], + b[2] - a[2], c[2] - a[2], d[2] - a[2]}; + + // det(M) = v1x(v2y*v3z - v2z*v3y) - v1y(v2x*v3z - v3x*v2z) + v1z(v2x*v3y - v3x*v2y) T determinant = (volumeMatrix[k_00] * (volumeMatrix[k_11] * volumeMatrix[k_22] - volumeMatrix[k_12] * volumeMatrix[k_21])) - - (volumeMatrix[k_01] * (volumeMatrix[k_10] * volumeMatrix[k_22] - volumeMatrix[k_12] * volumeMatrix[k_20])) + - (volumeMatrix[k_02] * (volumeMatrix[k_10] * volumeMatrix[k_21] - volumeMatrix[k_11] * volumeMatrix[k_20])); - - tetInfo[4 * iter + 0] = (determinant / 6.0f); - tetInfo[4 * iter + 1] = 0.25 * (ax + bx + cx + dx); - tetInfo[4 * iter + 2] = 0.25 * (ay + by + cy + dy); - tetInfo[4 * iter + 3] = 0.25 * (az + bz + cz + dz); + (volumeMatrix[k_10] * (volumeMatrix[k_01] * volumeMatrix[k_22] - volumeMatrix[k_02] * volumeMatrix[k_21])) + + (volumeMatrix[k_20] * (volumeMatrix[k_01] * volumeMatrix[k_12] - volumeMatrix[k_02] * volumeMatrix[k_11])); + // clang-format on + + tetInfo[4 * iter + 0] = std::abs((determinant / 6.0f)); // The final volume of the tetrahedron + // Encode the centroid of the sub-tetrahedron + tetInfo[4 * iter + 1] = 0.25 * (a[0] + b[0] + c[0] + d[0]); + tetInfo[4 * iter + 2] = 0.25 * (a[1] + b[1] + c[1] + d[1]); + tetInfo[4 * iter + 3] = 0.25 * (a[2] + b[2] + c[2] + d[2]); } } +using TriStore = AbstractDataStore; +using VertsStore = AbstractDataStore; + +/** + * @brief This will extract the 3 vertices from a given triangle face of a triangle geometry. This is MUCH faster + * than calling the function in the Triangle Geometry because of the dynamic_cast<> that goes on in that function. + */ +inline std::array GetFaceCoordinates(usize triangleId, const VertsStore& verts, const TriStore& triangleList) +{ + usize v0Idx = triangleList[triangleId * 3]; + usize v1Idx = triangleList[triangleId * 3 + 1]; + usize v2Idx = triangleList[triangleId * 3 + 2]; + return {Point3Df{verts[v0Idx * 3], verts[v0Idx * 3 + 1], verts[v0Idx * 3 + 2]}, Point3Df{verts[v1Idx * 3], verts[v1Idx * 3 + 1], verts[v1Idx * 3 + 2]}, + Point3Df{verts[v2Idx * 3], verts[v2Idx * 3 + 1], verts[v2Idx * 3 + 2]}}; +} + } // namespace // ----------------------------------------------------------------------------- void ComputeTriangleGeomShapes::findMoments() { using MeshIndexType = IGeometry::MeshIndexType; - using SharedVertexListType = IGeometry::SharedVertexList; + // using SharedVertexListType = IGeometry::SharedVertexList; - const auto& triangles = m_DataStructure.getDataRefAs(m_InputValues->TriangleGeometryPath); - const SharedVertexListType& vertPtr = triangles.getVerticesRef(); - const Int32Array& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath); - const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); - const Float32Array& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); + const auto& triangleGeom = m_DataStructure.getDataRefAs(m_InputValues->TriangleGeometryPath); + const TriStore& triangleList = triangleGeom.getFacesRef().getDataStoreRef(); + const VertsStore& verts = triangleGeom.getVerticesRef().getDataStoreRef(); + const auto& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath); + const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); + const auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); + // Calculated Arrays auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); usize numFaces = faceLabels.getNumberOfTuples(); @@ -121,30 +146,45 @@ void ComputeTriangleGeomShapes::findMoments() std::array centroid = {0.0F, 0.0F, 0.0F}; std::array tetInfo = {0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F}; - std::array vertIds = {0, 0, 0}; + // std::array vertIds = {0, 0, 0}; + // Loop over all triangle faces for(usize i = 0; i < numFaces; i++) { - triangles.getFacePointIds(i, vertIds); + // Get the indices to the 3 vertices of the triangle + // triangleGeom.getFacePointIds(i, vertIds); // This line is SLOW!!! + + std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); + + // Loop over each FaceLabel, which is basically the pair of "Feature Id" for the triangle for(usize j = 0; j < 2; j++) { + // Flips the winding, Assuming the triangle winding was correct in the first place. if(j == 1) { - std::swap(vertIds[2], vertIds[1]); + std::swap(vertCoords[2], vertCoords[1]); } int32 gnum = faceLabels[2 * i + j]; if(gnum > 0) { + // Get the centroid for the feature centroid[0] = centroids[3 * gnum + 0]; centroid[1] = centroids[3 * gnum + 1]; centroid[2] = centroids[3 * gnum + 2]; - FindTetrahedronInfo(vertIds, vertPtr, centroid, tetInfo); + // Generate the volumes for the 8 sub-tetrahedrons + FindTetrahedronInfo(vertCoords, centroid, tetInfo); + + // Loop over those 8 sub-tetrahedrons for(usize iter = 0; iter < 8; iter++) { + + // Compute the distance between the feature's centroid and the centroid of + // the current sub-tetrahedron. float64 xdist = (tetInfo[4 * iter + 1] - centroids[gnum * 3 + 0]); float64 ydist = (tetInfo[4 * iter + 2] - centroids[gnum * 3 + 1]); float64 zdist = (tetInfo[4 * iter + 3] - centroids[gnum * 3 + 2]); + // Compute the second order moments auto xx = static_cast((ydist * ydist) + (zdist * zdist)); auto yy = static_cast((xdist * xdist) + (zdist * zdist)); auto zz = static_cast((xdist * xdist) + (ydist * ydist)); @@ -152,6 +192,7 @@ void ComputeTriangleGeomShapes::findMoments() auto yz = static_cast(ydist * zdist); auto xz = static_cast(xdist * zdist); + // Sum the volume contributions into the "Feature Moments" array m_FeatureMoments[gnum * 6 + 0] = m_FeatureMoments[gnum * 6 + 0] + (xx * tetInfo[4 * iter + 0]); m_FeatureMoments[gnum * 6 + 1] = m_FeatureMoments[gnum * 6 + 1] + (yy * tetInfo[4 * iter + 0]); m_FeatureMoments[gnum * 6 + 2] = m_FeatureMoments[gnum * 6 + 2] + (zz * tetInfo[4 * iter + 0]); @@ -166,7 +207,7 @@ void ComputeTriangleGeomShapes::findMoments() const float64 k_Sphere = (2000.0 * M_PI * M_PI) / 9.0; for(usize i = 1; i < numFeatures; i++) { - float64 vol5 = pow(volumes[i], 5.0); + m_FeatureMoments[i * 6 + 3] = -m_FeatureMoments[i * 6 + 3]; m_FeatureMoments[i * 6 + 4] = -m_FeatureMoments[i * 6 + 4]; m_FeatureMoments[i * 6 + 5] = -m_FeatureMoments[i * 6 + 5]; @@ -177,6 +218,7 @@ void ComputeTriangleGeomShapes::findMoments() auto u011 = static_cast(-m_FeatureMoments[i * 6 + 4]); auto u101 = static_cast(-m_FeatureMoments[i * 6 + 5]); auto o3 = static_cast((u200 * u020 * u002) + (2.0f * u110 * u101 * u011) - (u200 * u011 * u011) - (u020 * u101 * u101) - (u002 * u110 * u110)); + float64 vol5 = pow(volumes[i], 5.0); float64 omega3 = vol5 / o3; omega3 = omega3 / k_Sphere; if(omega3 > 1) From b4b262ce9fbf49c4b95b58b73b4b487dc813fd25 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Fri, 6 Dec 2024 23:12:16 -0500 Subject: [PATCH 02/12] Compute Shape Values working. Axis and EulerAxis are not verified yet. Signed-off-by: Michael Jackson --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 321 ++++++------------ .../Algorithms/ComputeTriangleGeomShapes.hpp | 1 + src/simplnx/Utilities/GeometryUtilities.hpp | 47 ++- 3 files changed, 149 insertions(+), 220 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index a155ec6524..082b9a8a07 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -12,99 +12,6 @@ using namespace nx::core; namespace { -constexpr float64 k_ScaleFactor = 1.0; -constexpr usize k_00 = 0; -constexpr usize k_01 = 1; -constexpr usize k_02 = 2; - -constexpr usize k_10 = 3; -constexpr usize k_11 = 4; -constexpr usize k_12 = 5; - -constexpr usize k_20 = 6; -constexpr usize k_21 = 7; -constexpr usize k_22 = 8; - -// ----------------------------------------------------------------------------- -/** - * @brief This function will calculate the volume and centroid for 8 sub tetrahedrons that are generated - * from the triangle given by the "verts" and a centroid point give by "centroid" - * @tparam T - * @param verts The 3 Vertices of the triangle - * @param centroid Another point in space that represents the centroid of the feature that the triangle belongs to - * @param tetInfo The array of TetInfo data which is the volume and centroid of each sub-tetrahedron that is formed. - */ -template -void FindTetrahedronInfo(const std::array& verts, const std::array& centroid, std::array& tetInfo) -{ - // clang-format off - std::array coords = {verts[0][0], verts[0][1], verts[0][2], // V0 - verts[1][0], verts[1][1], verts[1][2], // V1 - verts[2][0], verts[2][1], verts[2][2], // V2 - - centroid[0], centroid[1], centroid[2], // Centroid Vertex - - 0.5 * (verts[0][0] + centroid[0]), 0.5 * (verts[0][1] + centroid[1]), 0.5 * (verts[0][2] + centroid[2]), // Average of V0 and Centroid - 0.5 * (verts[1][0] + centroid[0]), 0.5 * (verts[1][1] + centroid[1]), 0.5 * (verts[1][2] + centroid[2]), // Average of V1 and Centroid - 0.5 * (verts[2][0] + centroid[0]), 0.5 * (verts[2][1] + centroid[1]), 0.5 * (verts[2][2] + centroid[2]), // Average of V2 and Centroid - - 0.5 * (verts[0][0] + verts[1][0]), 0.5 * (verts[0][1] + verts[1][1]), 0.5 * (verts[0][2] + verts[1][2]), // Average of V0 & V1 Distance - 0.5 * (verts[1][0] + verts[2][0]), 0.5 * (verts[1][1] + verts[2][1]), 0.5 * (verts[1][2] + verts[2][2]), // Average of V1 & V2 Distance - 0.5 * (verts[2][0] + verts[0][0]), 0.5 * (verts[2][1] + verts[0][1]), 0.5 * (verts[2][2] + verts[0][2]) // Average of V2 & V0 Distance - }; - - // Create 8 Tetrahedrons - std::array tets = { - 4, 5, 6, 3, - 0, 7, 9, 4, - 1, 8, 7, 5, - 2, 9, 8, 6, - 7, 5, 6, 4, - 6, 9, 7, 4, - 6, 5, 7, 8, - 7, 9, 6, 8 - }; - // clang-format on - - using Point3DType = Point3D; - // Loop over each Tetrahedron - for(usize iter = 0; iter < 8; iter++) - { - usize i0 = 4 * iter + 0; - usize i1 = 4 * iter + 1; - usize i2 = 4 * iter + 2; - usize i3 = 4 * iter + 3; - // Create the 4 Vertices of the sub-tetrahedron - Point3DType a(coords[3 * tets[i0] + 0], coords[3 * tets[i0] + 1], coords[3 * tets[i0] + 2]); - Point3DType b(coords[3 * tets[i1] + 0], coords[3 * tets[i1] + 1], coords[3 * tets[i1] + 2]); - Point3DType c(coords[3 * tets[i2] + 0], coords[3 * tets[i2] + 1], coords[3 * tets[i2] + 2]); - Point3DType d(coords[3 * tets[i3] + 0], coords[3 * tets[i3] + 1], coords[3 * tets[i3] + 2]); - - // The volume of the tetrahedron can be calculated through V= 1/6 |M| (where M is the volume matrix) - // v1 = b-a, v2 = c-a, v3 = d-a - // | v1x v2x v3x - // M = | v1y v2y v3y - // | v1z v2z v3z - - // clang-format off - std::array volumeMatrix = {b[0] - a[0], c[0] - a[0], d[0] - a[0], - b[1] - a[1], c[1] - a[1], d[1] - a[1], - b[2] - a[2], c[2] - a[2], d[2] - a[2]}; - - // det(M) = v1x(v2y*v3z - v2z*v3y) - v1y(v2x*v3z - v3x*v2z) + v1z(v2x*v3y - v3x*v2y) - T determinant = (volumeMatrix[k_00] * (volumeMatrix[k_11] * volumeMatrix[k_22] - volumeMatrix[k_12] * volumeMatrix[k_21])) - - (volumeMatrix[k_10] * (volumeMatrix[k_01] * volumeMatrix[k_22] - volumeMatrix[k_02] * volumeMatrix[k_21])) + - (volumeMatrix[k_20] * (volumeMatrix[k_01] * volumeMatrix[k_12] - volumeMatrix[k_02] * volumeMatrix[k_11])); - // clang-format on - - tetInfo[4 * iter + 0] = std::abs((determinant / 6.0f)); // The final volume of the tetrahedron - // Encode the centroid of the sub-tetrahedron - tetInfo[4 * iter + 1] = 0.25 * (a[0] + b[0] + c[0] + d[0]); - tetInfo[4 * iter + 2] = 0.25 * (a[1] + b[1] + c[1] + d[1]); - tetInfo[4 * iter + 3] = 0.25 * (a[2] + b[2] + c[2] + d[2]); - } -} - using TriStore = AbstractDataStore; using VertsStore = AbstractDataStore; @@ -123,19 +30,45 @@ inline std::array GetFaceCoordinates(usize triangleId, co } // namespace +// ----------------------------------------------------------------------------- +ComputeTriangleGeomShapes::ComputeTriangleGeomShapes(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + ComputeTriangleGeomShapesInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeTriangleGeomShapes::~ComputeTriangleGeomShapes() noexcept = default; + +// ----------------------------------------------------------------------------- +const std::atomic_bool& ComputeTriangleGeomShapes::getCancel() +{ + return m_ShouldCancel; +} + +// ----------------------------------------------------------------------------- +Result<> ComputeTriangleGeomShapes::operator()() +{ + findMoments(); + findAxes(); + findAxisEulers(); + + return {}; +} + // ----------------------------------------------------------------------------- void ComputeTriangleGeomShapes::findMoments() { using MeshIndexType = IGeometry::MeshIndexType; - // using SharedVertexListType = IGeometry::SharedVertexList; - const auto& triangleGeom = m_DataStructure.getDataRefAs(m_InputValues->TriangleGeometryPath); const TriStore& triangleList = triangleGeom.getFacesRef().getDataStoreRef(); const VertsStore& verts = triangleGeom.getVerticesRef().getDataStoreRef(); const auto& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath); const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); - const auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); // Calculated Arrays auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); @@ -143,107 +76,97 @@ void ComputeTriangleGeomShapes::findMoments() usize numFeatures = centroids.getNumberOfTuples(); m_FeatureMoments.resize(numFeatures * 6, 0.0); - std::array centroid = {0.0F, 0.0F, 0.0F}; - std::array tetInfo = {0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, - 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F}; - // std::array vertIds = {0, 0, 0}; + nx::core::Point3Df centroid = {0.0F, 0.0F, 0.0F}; - // Loop over all triangle faces - for(usize i = 0; i < numFaces; i++) - { - // Get the indices to the 3 vertices of the triangle - // triangleGeom.getFacePointIds(i, vertIds); // This line is SLOW!!! + using Matrix3x3 = Eigen::Matrix; + // Theoretical perfect Sphere value of Omega-3. Each calculated Omega-3 + // will be normalized using this value; + const float64 k_Sphere = (2000.0 * M_PI * M_PI) / 9.0; - std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); + // define the canonical C matrix + double aa = 1.0 / 60.0; + double bb = aa / 2.0; + // clang-format off + Matrix3x3 C; + C << aa, bb, bb, bb, aa, bb, bb, bb, aa; + + // and the identity matrix + aa = 1.0; + bb = 0.0; + Matrix3x3 ID; + ID << aa, bb, bb, bb, aa, bb, bb, bb, aa; + + // The C-Prime matrix + Matrix3x3 CC; + CC << -0.50000000, 0.50000000, 0.50000000, + 0.50000000, -0.50000000, 0.50000000, + 0.50000000, 0.50000000, -0.50000000; + // clang-format on - // Loop over each FaceLabel, which is basically the pair of "Feature Id" for the triangle - for(usize j = 0; j < 2; j++) + // Loop over each "Feature" which is the number of tuples in the "Centroids" array + // We could parallelize over the features? + for(usize featureId = 1; featureId < numFeatures; featureId++) + { + if(m_ShouldCancel) { - // Flips the winding, Assuming the triangle winding was correct in the first place. - if(j == 1) + return; + } + double Vol = 0.0; + // define the accumulator arrays + Matrix3x3 Cacc; + Cacc << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0; + + // Get the centroid for the feature + centroid[0] = centroids[3 * featureId + 0]; + centroid[1] = centroids[3 * featureId + 1]; + centroid[2] = centroids[3 * featureId + 2]; + + // for each triangle we need the transformation matrix A defined by the three points as columns + // Loop over all triangle faces + int32_t tCount = 0; + for(usize i = 0; i < numFaces; i++) + { + if(faceLabels[2 * i] != featureId && faceLabels[2 * i + 1] != featureId) { - std::swap(vertCoords[2], vertCoords[1]); + continue; } - int32 gnum = faceLabels[2 * i + j]; - if(gnum > 0) - { - // Get the centroid for the feature - centroid[0] = centroids[3 * gnum + 0]; - centroid[1] = centroids[3 * gnum + 1]; - centroid[2] = centroids[3 * gnum + 2]; - // Generate the volumes for the 8 sub-tetrahedrons - FindTetrahedronInfo(vertCoords, centroid, tetInfo); - - // Loop over those 8 sub-tetrahedrons - for(usize iter = 0; iter < 8; iter++) - { + tCount++; + usize compIndex = (faceLabels[2 * i] == featureId ? 0 : 1); + std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); - // Compute the distance between the feature's centroid and the centroid of - // the current sub-tetrahedron. - float64 xdist = (tetInfo[4 * iter + 1] - centroids[gnum * 3 + 0]); - float64 ydist = (tetInfo[4 * iter + 2] - centroids[gnum * 3 + 1]); - float64 zdist = (tetInfo[4 * iter + 3] - centroids[gnum * 3 + 2]); - - // Compute the second order moments - auto xx = static_cast((ydist * ydist) + (zdist * zdist)); - auto yy = static_cast((xdist * xdist) + (zdist * zdist)); - auto zz = static_cast((xdist * xdist) + (ydist * ydist)); - auto xy = static_cast(xdist * ydist); - auto yz = static_cast(ydist * zdist); - auto xz = static_cast(xdist * zdist); - - // Sum the volume contributions into the "Feature Moments" array - m_FeatureMoments[gnum * 6 + 0] = m_FeatureMoments[gnum * 6 + 0] + (xx * tetInfo[4 * iter + 0]); - m_FeatureMoments[gnum * 6 + 1] = m_FeatureMoments[gnum * 6 + 1] + (yy * tetInfo[4 * iter + 0]); - m_FeatureMoments[gnum * 6 + 2] = m_FeatureMoments[gnum * 6 + 2] + (zz * tetInfo[4 * iter + 0]); - m_FeatureMoments[gnum * 6 + 3] = m_FeatureMoments[gnum * 6 + 3] + (xy * tetInfo[4 * iter + 0]); - m_FeatureMoments[gnum * 6 + 4] = m_FeatureMoments[gnum * 6 + 4] + (yz * tetInfo[4 * iter + 0]); - m_FeatureMoments[gnum * 6 + 5] = m_FeatureMoments[gnum * 6 + 5] + (xz * tetInfo[4 * iter + 0]); - } - } - } - } + const nx::core::Point3Df& a = vertCoords[0] - centroid; + const nx::core::Point3Df& b = (compIndex == 0 ? vertCoords[1] : vertCoords[2]) - centroid; + const nx::core::Point3Df& c = (compIndex == 0 ? vertCoords[2] : vertCoords[1]) - centroid; - const float64 k_Sphere = (2000.0 * M_PI * M_PI) / 9.0; - for(usize i = 1; i < numFeatures; i++) - { + Matrix3x3 A; + A << a[0], b[0], c[0], a[1], b[1], c[1], a[2], b[2], c[2]; - m_FeatureMoments[i * 6 + 3] = -m_FeatureMoments[i * 6 + 3]; - m_FeatureMoments[i * 6 + 4] = -m_FeatureMoments[i * 6 + 4]; - m_FeatureMoments[i * 6 + 5] = -m_FeatureMoments[i * 6 + 5]; - auto u200 = static_cast((m_FeatureMoments[i * 6 + 1] + m_FeatureMoments[i * 6 + 2] - m_FeatureMoments[i * 6 + 0]) / 2.0f); - auto u020 = static_cast((m_FeatureMoments[i * 6 + 0] + m_FeatureMoments[i * 6 + 2] - m_FeatureMoments[i * 6 + 1]) / 2.0f); - auto u002 = static_cast((m_FeatureMoments[i * 6 + 0] + m_FeatureMoments[i * 6 + 1] - m_FeatureMoments[i * 6 + 2]) / 2.0f); - auto u110 = static_cast(-m_FeatureMoments[i * 6 + 3]); - auto u011 = static_cast(-m_FeatureMoments[i * 6 + 4]); - auto u101 = static_cast(-m_FeatureMoments[i * 6 + 5]); - auto o3 = static_cast((u200 * u020 * u002) + (2.0f * u110 * u101 * u011) - (u200 * u011 * u011) - (u020 * u101 * u101) - (u002 * u110 * u110)); - float64 vol5 = pow(volumes[i], 5.0); - float64 omega3 = vol5 / o3; - omega3 = omega3 / k_Sphere; - if(omega3 > 1) - { - omega3 = 1.0; - } - if(vol5 == 0.0) - { - omega3 = 0.0; + float64 dA = A.determinant(); + + Cacc = (Cacc + dA * (A * (C * (A.transpose())))).eval(); + Vol += (dA / 6.0f); // accumulate the volumes } - omega3S[i] = static_cast(omega3); + + Cacc = (Cacc / Vol).eval(); + Matrix3x3 Cinertia = ID * Cacc.trace() - Cacc; + // extract the moments from the inertia tensor + Eigen::Vector3d e(Cinertia(0, 0), Cinertia(1, 1), Cinertia(2, 2)); + auto sols = CC * e; + omega3S[featureId] = static_cast(((Vol * Vol) / sols.prod()) / k_Sphere); + + m_FeatureMoments[featureId * 6 + 0] = sols[0]; + m_FeatureMoments[featureId * 6 + 1] = sols[1]; + m_FeatureMoments[featureId * 6 + 2] = sols[2]; + m_FeatureMoments[featureId * 6 + 3] = -Cinertia(0, 1); + m_FeatureMoments[featureId * 6 + 4] = -Cinertia(0, 2); + m_FeatureMoments[featureId * 6 + 5] = -Cinertia(1, 2); } } // ----------------------------------------------------------------------------- void ComputeTriangleGeomShapes::findAxes() { - // float64 I1 = 0.0, I2 = 0.0, I3 = 0.0; - // float64 a = 0.0, b = 0.0, c = 0.0, d = 0.0, f = 0.0, g = 0.0, h = 0.0; - // float64 rsquare = 0.0, r = 0.0, theta = 0.0; - // float64 A = 0.0, B = 0.0, C = 0.0; - // float64 r1 = 0.0, r2 = 0.0, r3 = 0.0; - // float32 bovera = 0.0f, covera = 0.0f; - // float64 value = 0.0; - + constexpr float64 k_ScaleFactor = 1.0; const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); @@ -255,6 +178,10 @@ void ComputeTriangleGeomShapes::findAxes() for(usize i = 1; i < numFeatures; i++) { + if(m_ShouldCancel) + { + return; + } float64 ixx = m_FeatureMoments[i * 6 + 0]; float64 iyy = m_FeatureMoments[i * 6 + 1]; float64 izz = m_FeatureMoments[i * 6 + 2]; @@ -344,6 +271,10 @@ void ComputeTriangleGeomShapes::findAxisEulers() for(usize i = 1; i < numFeatures; i++) { + if(m_ShouldCancel) + { + return; + } float64 ixx = m_FeatureMoments[i * 6 + 0]; float64 iyy = m_FeatureMoments[i * 6 + 1]; float64 izz = m_FeatureMoments[i * 6 + 2]; @@ -470,31 +401,3 @@ void ComputeTriangleGeomShapes::findAxisEulers() } } -// ----------------------------------------------------------------------------- -ComputeTriangleGeomShapes::ComputeTriangleGeomShapes(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, - ComputeTriangleGeomShapesInputValues* inputValues) -: m_DataStructure(dataStructure) -, m_InputValues(inputValues) -, m_ShouldCancel(shouldCancel) -, m_MessageHandler(mesgHandler) -{ -} - -// ----------------------------------------------------------------------------- -ComputeTriangleGeomShapes::~ComputeTriangleGeomShapes() noexcept = default; - -// ----------------------------------------------------------------------------- -const std::atomic_bool& ComputeTriangleGeomShapes::getCancel() -{ - return m_ShouldCancel; -} - -// ----------------------------------------------------------------------------- -Result<> ComputeTriangleGeomShapes::operator()() -{ - findMoments(); - findAxes(); - findAxisEulers(); - - return {}; -} diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp index 839e83cb53..d48bd1f139 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp @@ -48,6 +48,7 @@ class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapes std::vector m_FeatureMoments; std::vector m_FeatureEigenVals; + void findMoments(); void findAxes(); void findAxisEulers(); diff --git a/src/simplnx/Utilities/GeometryUtilities.hpp b/src/simplnx/Utilities/GeometryUtilities.hpp index 29fea11bcf..1b44e8e183 100644 --- a/src/simplnx/Utilities/GeometryUtilities.hpp +++ b/src/simplnx/Utilities/GeometryUtilities.hpp @@ -64,17 +64,42 @@ SIMPLNX_EXPORT Result CalculateNodeBasedPartitionSchemeOrigin(const I */ SIMPLNX_EXPORT Result CalculatePartitionLengthsOfBoundingBox(const BoundingBox3Df& boundingBox, const SizeVec3& numberOfPartitionsPerAxis); -///** -// * @brief Constructs a new BoundingBox3D defined by the array of position values. -// * The format is min X, min Y, min Z, max X, max Y, max Z. -// * @param arr -// */ -// template >::value>> -// explicit BoundingBox(nonstd::span arr) -//: m_Lower(Point3D(arr[0], arr[1], arr[2])) -//, m_Upper(Point3D(arr[3], arr[4], arr[5])) -//{ -//} +/** + * @brief This function will find the volume of a tetrahedron + * @tparam T + * @param verts + * @param v4 + * @return + */ +template +T ComputeTetrahedronVolume(const std::array& verts, const std::array& v4) +{ + using Point3DType = Point3D; + + // Create the 4 Vertices of the tetrahedron + const Point3DType& a = verts[0]; + const Point3DType& b = verts[1]; + const Point3DType& c = verts[2]; + const Point3DType& d = v4; + + // The volume of the tetrahedron can be calculated through V= 1/6 |M| (where M is the volume matrix) + // v1 = b-a, v2 = c-a, v3 = d-a + // | v1x v2x v3x | + // M = | v1y v2y v3y | + // | v1z v2z v3z | + // clang-format off + std::array volumeMatrix = {b[0] - a[0], c[0] - a[0], d[0] - a[0], + b[1] - a[1], c[1] - a[1], d[1] - a[1], + b[2] - a[2], c[2] - a[2], d[2] - a[2]}; + + // det(M) = v1x(v2y*v3z - v2z*v3y) - v1y(v2x*v3z - v3x*v2z) + v1z(v2x*v3y - v3x*v2y) + T determinant = (volumeMatrix[0] * (volumeMatrix[4] * volumeMatrix[8] - volumeMatrix[5] * volumeMatrix[7])) - + (volumeMatrix[3] * (volumeMatrix[1] * volumeMatrix[8] - volumeMatrix[2] * volumeMatrix[7])) + + (volumeMatrix[6] * (volumeMatrix[1] * volumeMatrix[5] - volumeMatrix[2] * volumeMatrix[4])); + // clang-format on + + return std::abs((determinant / 6.0f)); // The final volume of the tetrahedron +} /** * @brief Removes duplicate nodes to ensure the vertex list is unique From 71986c9afc6f0c6b17cb24743c0e92a9b9a0d613 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Sat, 7 Dec 2024 01:03:47 -0500 Subject: [PATCH 03/12] Implement Unit test for Omega-3 Values only. All values are validated against the voxel generated results and the theoretical values. Theoretical values taken from paper from De Graef, Simmons, et. al. Signed-off-by: Michael Jackson --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 1 - .../OrientationAnalysis/test/CMakeLists.txt | 2 +- .../test/ComputeTriangleGeomShapesTest.cpp | 137 +++++++++++------- .../test/MergeTwinsTest.cpp | 2 - 4 files changed, 83 insertions(+), 59 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 082b9a8a07..7486104452 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -400,4 +400,3 @@ void ComputeTriangleGeomShapes::findAxisEulers() axisEulerAngles[3 * i + 2] = euler[2]; } } - diff --git a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt index aa13c0a38d..8729c9ee5d 100644 --- a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt +++ b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt @@ -31,7 +31,7 @@ set(${PLUGIN_NAME}UnitTest_SRCS ComputeSchmidsTest.cpp ComputeShapesFilterTest.cpp ComputeSlipTransmissionMetricsTest.cpp - # ComputeTriangleGeomShapesTest.cpp + ComputeTriangleGeomShapesTest.cpp ConvertHexGridToSquareGridTest.cpp ConvertOrientationsTest.cpp ConvertQuaternionTest.cpp diff --git a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp index 473efe0c44..776e6659df 100644 --- a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp @@ -8,15 +8,22 @@ #include +#include +#include +#include + +namespace fs = std::filesystem; + using namespace nx::core; using namespace nx::core::UnitTest; +#define SIMPLNX_WRITE_TEST_OUTPUT + namespace ComputeTriangleGeomShapesFilterTest { -const std::string k_TriangleGeometryName = "TriangleDataContainer"; -const std::string k_FaceLabelsName = "FaceLabels"; -const std::string k_FaceFeatureName = "FaceFeatureData"; -const std::string k_FaceDataName = "FaceData"; +const std::string k_FaceLabelsName = "Face Labels"; +const std::string k_FaceFeatureName = "Face Feature Data"; +const std::string k_FaceDataName = "Face Data"; const std::string k_CentroidsArrayName = "Centroids"; const std::string k_VolumesArrayName = "Volumes"; @@ -25,11 +32,31 @@ const std::string k_AxisLengthsArrayName = "AxisLengths [NX Computed]"; const std::string k_AxisEulerAnglesArrayName = "AxisEulerAngles [NX Computed]"; const std::string k_AspectRatiosArrayName = "AspectRatios [NX Computed]"; -const DataPath k_GeometryPath = DataPath({k_TriangleGeometryName}); -const DataPath k_FaceFeatureAttributeMatrixPath = k_GeometryPath.createChildPath(k_FaceFeatureName); -const DataPath k_FaceLabelsPath = k_GeometryPath.createChildPath(k_FaceDataName).createChildPath(k_FaceLabelsName); -const DataPath k_FaceFeatureCentroidsPath = k_FaceFeatureAttributeMatrixPath.createChildPath(k_CentroidsArrayName); -const DataPath k_FaceFeatureVolumesPath = k_FaceFeatureAttributeMatrixPath.createChildPath(k_VolumesArrayName); +// const DataPath k_FaceFeatureAttributeMatrixPath = k_GeometryPath.createChildPath(k_FaceFeatureName); +// const DataPath k_FaceLabelsPath = k_GeometryPath.createChildPath(k_FaceDataName).createChildPath(k_FaceLabelsName); +// const DataPath k_FaceFeatureCentroidsPath = k_FaceFeatureAttributeMatrixPath.createChildPath(k_CentroidsArrayName); +// const DataPath k_FaceFeatureVolumesPath = k_FaceFeatureAttributeMatrixPath.createChildPath(k_VolumesArrayName); + +constexpr float32 k_Sphere_Omega_3 = 2193.245f; // 1.00000 + +constexpr float32 k_Tetrahedron_Omega_3 = 888.889f / k_Sphere_Omega_3; +constexpr float32 k_Cube_Omega_3 = 1728.0f / k_Sphere_Omega_3; // 0.787873675763538 +constexpr float32 k_Octahedron_Omega_3 = 1777.778f / k_Sphere_Omega_3; // 0.810569726592332 +constexpr float32 k_Dodecahedron_Omega_3 = 2096.873f / k_Sphere_Omega_3; // 0.956059628541271 +constexpr float32 k_Icosahedron_Omega_3 = 2122.033f / k_Sphere_Omega_3; // 0.967531215162921 +constexpr float32 k_TriangularPyramid_Omega_3 = 888.889f / k_Sphere_Omega_3; // 0.405284863296166 +constexpr float32 k_Rectangular_Prism_Omega_3 = 1728.0f / k_Sphere_Omega_3; // 0.787873675763538 +constexpr float32 k_Circular_Cylinder_Omega_3 = 1894.964f / k_Sphere_Omega_3; // 0.864000145902533 +constexpr float32 k_Ellipsoid_Omega_3 = 2193.245f / k_Sphere_Omega_3; // 1.000 + +const std::string k_TestFileDirName = "7_Triangle_Shapes_Files/stl_exemplar"; +std::vector k_TestFileNames = {"10_octahedron", "11_pyramid_elongated", "20_sphere", "21_ellipsoid", "40_rounded_cube", "41_rounded_cube_elongated", + "50_cube", "51_rectangular_prism", "60_Cylinder", "70_Dodecahedron", "80_Icosahedron", "90_tetrahedron"}; + +std::vector k_Theoretical_Omega_3_Values = { + k_Octahedron_Omega_3, k_Octahedron_Omega_3, 1.0f, k_Ellipsoid_Omega_3, 0.907957, 0.907945, k_Cube_Omega_3, k_Rectangular_Prism_Omega_3, k_Circular_Cylinder_Omega_3, k_Dodecahedron_Omega_3, + k_Icosahedron_Omega_3, k_Tetrahedron_Omega_3}; + } // namespace ComputeTriangleGeomShapesFilterTest using namespace ComputeTriangleGeomShapesFilterTest; @@ -37,54 +64,54 @@ using namespace ComputeTriangleGeomShapesFilterTest; TEST_CASE("OrientationAnalysis::ComputeTriangleGeomShapes", "[OrientationAnalysis][ComputeTriangleGeomShapes]") { Application::GetOrCreateInstance()->loadPlugins(unit_test::k_BuildDir.view(), true); - - const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "12_IN625_GBCD.tar.gz", "12_IN625_GBCD"); - - // Read Exemplar DREAM3D File Filter - auto exemplarFilePath = fs::path(fmt::format("{}/12_IN625_GBCD/12_IN625_GBCD.dream3d", unit_test::k_TestFilesDir)); - DataStructure dataStructure = LoadDataStructure(exemplarFilePath); - + const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "7_Triangle_Shapes_Files.tar.gz", "7_Triangle_Shapes_Files"); + usize index = 0; + for(const auto& fileName : k_TestFileNames) { - // Instantiate the filter and an Arguments Object - ComputeTriangleGeomShapesFilter filter; - Arguments args; - - // Create default Parameters for the filter. - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_TriGeometryDataPath_Key, std::make_any(k_GeometryPath)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_FaceLabelsArrayPath_Key, std::make_any(k_FaceLabelsPath)); - - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_FeatureAttributeMatrixPath_Key, std::make_any(k_FaceFeatureAttributeMatrixPath)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_CentroidsArrayPath_Key, std::make_any(k_FaceFeatureCentroidsPath)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_VolumesArrayPath_Key, std::make_any(k_FaceFeatureVolumesPath)); - // Output Vars - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_Omega3sArrayName_Key, std::make_any(k_Omega3SArrayName)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AxisLengthsArrayName_Key, std::make_any(k_AxisLengthsArrayName)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AxisEulerAnglesArrayName_Key, std::make_any(k_AxisEulerAnglesArrayName)); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AspectRatiosArrayName_Key, std::make_any(k_AspectRatiosArrayName)); - - // Preflight the filter and check result - auto preflightResult = filter.preflight(dataStructure, args); - SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); - - // Execute the filter and check the result - auto executeResult = filter.execute(dataStructure, args); - SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); - } - - std::vector outputArrayNames = {k_Omega3SArrayName, k_AxisLengthsArrayName, k_AxisEulerAnglesArrayName, k_AspectRatiosArrayName}; - std::vector exemplarArrayNames = {"Omega3s", "AxisLengths", "AxisEulerAngles", "AspectRatios"}; - for(usize i = 0; i < 4; i++) - { - const DataPath kExemplarArrayPath = k_FaceFeatureAttributeMatrixPath.createChildPath(exemplarArrayNames[i]); - const DataPath kNxArrayPath = k_FaceFeatureAttributeMatrixPath.createChildPath(outputArrayNames[i]); - - const auto& kExemplarsArray = dataStructure.getDataRefAs(kExemplarArrayPath); - const auto& kNxArray = dataStructure.getDataRefAs(kNxArrayPath); - - UnitTest::CompareDataArrays(kExemplarsArray, kNxArray); - } + // Read Exemplar DREAM3D File Filter + auto exemplarFilePath = fs::path(fmt::format("{}/{}/{}.dream3d", nx::core::unit_test::k_TestFilesDir, k_TestFileDirName, fileName)); + // std::cout << "Loading File: " << exemplarFilePath.string() << "\n"; + DataStructure dataStructure = UnitTest::LoadDataStructure(exemplarFilePath); + + DataPath geometryPath({fileName}); + DataPath faceDataPath = geometryPath.createChildPath(k_FaceDataName); + DataPath faceFeatureDataPath = geometryPath.createChildPath(k_FaceFeatureName); + + { + // Instantiate the filter and an Arguments Object + ComputeTriangleGeomShapesFilter filter; + Arguments args; + + // Create default Parameters for the filter. + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_TriGeometryDataPath_Key, std::make_any(geometryPath)); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_FaceLabelsArrayPath_Key, std::make_any(faceDataPath.createChildPath(k_FaceLabelsName))); + + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_FeatureAttributeMatrixPath_Key, std::make_any(faceFeatureDataPath)); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_CentroidsArrayPath_Key, std::make_any(faceFeatureDataPath.createChildPath(k_CentroidsArrayName))); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_VolumesArrayPath_Key, std::make_any(faceFeatureDataPath.createChildPath(k_VolumesArrayName))); + // Output Vars + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_Omega3sArrayName_Key, std::make_any(k_Omega3SArrayName)); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AxisLengthsArrayName_Key, std::make_any(k_AxisLengthsArrayName)); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AxisEulerAnglesArrayName_Key, std::make_any(k_AxisEulerAnglesArrayName)); + args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AspectRatiosArrayName_Key, std::make_any(k_AspectRatiosArrayName)); + + // Preflight the filter and check result + auto preflightResult = filter.preflight(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + + // Execute the filter and check the result + auto executeResult = filter.execute(dataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + } #ifdef SIMPLNX_WRITE_TEST_OUTPUT - WriteTestDataStructure(dataStructure, fs::path(fmt::format("{}/find_triangle_geom_shapes.dream3d", unit_test::k_BinaryTestOutputDir))); + WriteTestDataStructure(dataStructure, fs::path(fmt::format("{}/{}.dream3d", unit_test::k_BinaryTestOutputDir, fileName))); #endif + + const auto& computedOmega3 = dataStructure.getDataRefAs(faceFeatureDataPath.createChildPath(k_Omega3SArrayName)); + float32 diff = std::abs(computedOmega3[1] - k_Theoretical_Omega_3_Values[index]); + // std::cout << fileName << ": Omega-3 Computed: " << computedOmega3[1] << " Theoretical Value: " << k_Theoretical_Omega_3_Values[index] << " Diff: " << diff << "\n"; + REQUIRE(diff < 1.0E-4); + index++; + } } diff --git a/src/Plugins/OrientationAnalysis/test/MergeTwinsTest.cpp b/src/Plugins/OrientationAnalysis/test/MergeTwinsTest.cpp index fd0baea0b5..1b78b0bf7a 100644 --- a/src/Plugins/OrientationAnalysis/test/MergeTwinsTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/MergeTwinsTest.cpp @@ -23,8 +23,6 @@ TEST_CASE("Reconstruction::MergeTwinsFilter: Valid Execution", "[Reconstruction] const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "neighbor_orientation_correlation.tar.gz", "neighbor_orientation_correlation.dream3d"); - - Application::GetOrCreateInstance()->loadPlugins(unit_test::k_BuildDir.view(), true); auto* filterList = Application::Instance()->getFilterList(); // Make sure we can load the needed filters from the plugins From 87be91e6896226de7a19cfc3ac802030f6f88b28 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Fri, 20 Dec 2024 09:49:56 -0500 Subject: [PATCH 04/12] Use eigenvalues and eigenvector to find principle axes and primary axes --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 48 +++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 7486104452..5c8d985075 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -160,6 +160,54 @@ void ComputeTriangleGeomShapes::findMoments() m_FeatureMoments[featureId * 6 + 3] = -Cinertia(0, 1); m_FeatureMoments[featureId * 6 + 4] = -Cinertia(0, 2); m_FeatureMoments[featureId * 6 + 5] = -Cinertia(1, 2); + + /** + * This next section finds the principle axis via eigenvalues. + * Paper/Lecture Notes (Page 5): https://ocw.mit.edu/courses/16-07-dynamics-fall-2009/dd277ec654440f4c2b5b07d6c286c3fd_MIT16_07F09_Lec26.pdf + * Video Walkthrough [0:00-10:45]: https://www.youtube.com/watch?v=IEDniK9kmaw + * + * The main goal is to derive the eigenvalues from the moment of inertia tensor therein finding the eigenvectors, + * which are the angular velocity vectors. + * + * Code Keynotes for reviewers: + * - Hessenburg Decomposition is pre-processing to get an upper/right triangular matrix. + * - This significantly reduces the iterations needed for QR Decomposition + * - QR Factorization is different than QR Decomposition. + * - QR Decomposition expresses the product of an Orthogonal Matrix (Q) and an right/upper triangular matrix (R) as a singular matrix (A). + */ + // TODO: + // - Remove Copying between processing (inline integration) + // - Extract real part of complexes stored in eigenvalues/vectors + Eigen::HessenbergDecomposition hessDecomp(Cinertia); + Matrix3x3 hessenMatrixUpper = hessDecomp.matrixH(); + + Eigen::HouseholderQR hQR(hessenMatrixUpper); + Matrix3x3 hqrMatrix = hQR.matrixQR(); + + // Extract eigenvalues and eigenvectors + Eigen::EigenSolver eigenSolver(hqrMatrix); + + // The primary axis is the largest eigenvector + Eigen::EigenSolver::EigenvalueType eigenvalues = eigenSolver.eigenvalues(); + + // This is the angular velocity vector, each row represents an axial alignment (principle axis) + Eigen::EigenSolver::EigenvectorsType eigenvectors = eigenSolver.eigenvectors(); + + 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; } } From db939ddc6bc6dfaf0504fb33edf04f02d3564d52 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Fri, 20 Dec 2024 16:49:02 -0500 Subject: [PATCH 05/12] correct production of eigenvalues/vectors --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 53 +++++++++++++------ 1 file changed, 36 insertions(+), 17 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 5c8d985075..84bba9957f 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -175,39 +175,58 @@ void ComputeTriangleGeomShapes::findMoments() * - QR Factorization is different than QR Decomposition. * - QR Decomposition expresses the product of an Orthogonal Matrix (Q) and an right/upper triangular matrix (R) as a singular matrix (A). */ + + // Zero out small numbers in Cinertia: https://stackoverflow.com/a/54505281 + Cinertia = (0.000000001 < Cinertia.array().abs()).select(Cinertia, 0.0); + // TODO: // - Remove Copying between processing (inline integration) // - Extract real part of complexes stored in eigenvalues/vectors Eigen::HessenbergDecomposition hessDecomp(Cinertia); Matrix3x3 hessenMatrixUpper = hessDecomp.matrixH(); - Eigen::HouseholderQR hQR(hessenMatrixUpper); + // It is important to maintain the FullPivHouseholderQR implementation here, it is rank preserving and most numerically stable + Eigen::FullPivHouseholderQR hQR(hessenMatrixUpper); Matrix3x3 hqrMatrix = hQR.matrixQR(); // Extract eigenvalues and eigenvectors Eigen::EigenSolver eigenSolver(hqrMatrix); - // The primary axis is the largest eigenvector + // The primary axis is the largest eigenvalue Eigen::EigenSolver::EigenvalueType eigenvalues = eigenSolver.eigenvalues(); // This is the angular velocity vector, each row represents an axial alignment (principle axis) Eigen::EigenSolver::EigenvectorsType eigenvectors = eigenSolver.eigenvectors(); - 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; + /** + * Following section for debugging + */ + // 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; + // + // // Condition Number calculations from here: https://stackoverflow.com/questions/33575478/how-can-you-find-the-condition-number-in-eigen + // // Calculate Condition Number from JacobiSVD + // Eigen::JacobiSVD svd(hqrMatrix); + // double cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size() - 1); + // std::cout << "\n Condition Number from JacobiSVD: " << cond << std::endl; + // + // // Calculate Condition Number from Norms + // double condNumNorms = hqrMatrix.completeOrthogonalDecomposition().pseudoInverse().norm() * hqrMatrix.norm(); + // + // std::cout << "\n Condition Number from Norms: " << condNumNorms << std::endl; } } From fe2c3d33176fb7787dcfe5ffb9a562d6f6992a50 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Tue, 31 Dec 2024 15:48:47 -0500 Subject: [PATCH 06/12] Correct axis length production, spacial requirement reduction, code consolidation, further comment documentation --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 497 ++++++++---------- .../Algorithms/ComputeTriangleGeomShapes.hpp | 8 - 2 files changed, 205 insertions(+), 300 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 84bba9957f..65ae8e05f1 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -1,5 +1,6 @@ #include "ComputeTriangleGeomShapes.hpp" +#include "simplnx/Common/Constants.hpp" #include "simplnx/DataStructure/DataArray.hpp" #include "simplnx/DataStructure/DataGroup.hpp" #include "simplnx/DataStructure/Geometry/IGeometry.hpp" @@ -15,6 +16,9 @@ namespace using TriStore = AbstractDataStore; using VertsStore = AbstractDataStore; +constexpr double k_Multiplier = 1.0 / (4.0 * Constants::k_PiD); +constexpr float64 k_ScaleFactor = 1.0; + /** * @brief This will extract the 3 vertices from a given triangle face of a triangle geometry. This is MUCH faster * than calling the function in the Triangle Geometry because of the dynamic_cast<> that goes on in that function. @@ -28,6 +32,81 @@ inline std::array GetFaceCoordinates(usize triangleId, co Point3Df{verts[v2Idx * 3], verts[v2Idx * 3 + 1], verts[v2Idx * 3 + 2]}}; } +/** + * @brief Sorts the 3 values + * @param a First Value + * @param b Second Value + * @param c Third Value + * @return The indices in their sorted order + */ +template +std::array TripletSort(T a, T b, T c, bool lowToHigh) +{ + constexpr size_t A = 0; + constexpr size_t B = 1; + constexpr size_t C = 2; + std::array idx = {0, 1, 2}; + if(a > b && a > c) + { + // sorted[2] = a; + if(b > c) + { + // sorted[1] = b; + // sorted[0] = c; + idx = {C, B, A}; + } + else + { + // sorted[1] = c; + // sorted[0] = b; + idx = {B, C, A}; + } + } + else if(b > a && b > c) + { + // sorted[2] = b; + if(a > c) + { + // sorted[1] = a; + // sorted[0] = c; + idx = {C, A, B}; + } + else + { + // sorted[1] = c; + // sorted[0] = a; + idx = {A, C, B}; + } + } + else if(a > b) + { + // sorted[1] = a; + // sorted[0] = b; + // sorted[2] = c; + idx = {B, A, C}; + } + else if(a >= c && b >= c) + { + // sorted[0] = c; + // sorted[1] = a; + // sorted[2] = b; + idx = {C, A, B}; + } + else + { + // sorted[0] = a; + // sorted[1] = b; + // sorted[2] = c; + idx = {A, B, C}; + } + + if(!lowToHigh) + { + std::swap(idx[0], idx[2]); + } + return idx; +} + } // namespace // ----------------------------------------------------------------------------- @@ -51,16 +130,6 @@ const std::atomic_bool& ComputeTriangleGeomShapes::getCancel() // ----------------------------------------------------------------------------- Result<> ComputeTriangleGeomShapes::operator()() -{ - findMoments(); - findAxes(); - findAxisEulers(); - - return {}; -} - -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findMoments() { using MeshIndexType = IGeometry::MeshIndexType; const auto& triangleGeom = m_DataStructure.getDataRefAs(m_InputValues->TriangleGeometryPath); @@ -69,16 +138,21 @@ void ComputeTriangleGeomShapes::findMoments() const auto& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath); const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); + // Calculated Arrays auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); + auto& axisEulerAngles = m_DataStructure.getDataRefAs(m_InputValues->AxisEulerAnglesArrayPath); + auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); + auto& aspectRatios = m_DataStructure.getDataRefAs(m_InputValues->AspectRatiosArrayPath); + + using Matrix3x3 = Eigen::Matrix; + Matrix3x3 Cinertia; usize numFaces = faceLabels.getNumberOfTuples(); usize numFeatures = centroids.getNumberOfTuples(); - m_FeatureMoments.resize(numFeatures * 6, 0.0); nx::core::Point3Df centroid = {0.0F, 0.0F, 0.0F}; - using Matrix3x3 = Eigen::Matrix; // Theoretical perfect Sphere value of Omega-3. Each calculated Omega-3 // will be normalized using this value; const float64 k_Sphere = (2000.0 * M_PI * M_PI) / 9.0; @@ -107,59 +181,57 @@ void ComputeTriangleGeomShapes::findMoments() // We could parallelize over the features? for(usize featureId = 1; featureId < numFeatures; featureId++) { - if(m_ShouldCancel) - { - return; - } - double Vol = 0.0; - // define the accumulator arrays - Matrix3x3 Cacc; - Cacc << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0; - - // Get the centroid for the feature - centroid[0] = centroids[3 * featureId + 0]; - centroid[1] = centroids[3 * featureId + 1]; - centroid[2] = centroids[3 * featureId + 2]; - - // for each triangle we need the transformation matrix A defined by the three points as columns - // Loop over all triangle faces - int32_t tCount = 0; - for(usize i = 0; i < numFaces; i++) + /** + * The following section calculates moment of inertia tensor (Cinertia) and omega3s + */ { - if(faceLabels[2 * i] != featureId && faceLabels[2 * i + 1] != featureId) + if(m_ShouldCancel) { - continue; + return {}; } - tCount++; - usize compIndex = (faceLabels[2 * i] == featureId ? 0 : 1); - std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); + double Vol = 0.0; + // define the accumulator arrays + Matrix3x3 Cacc; + Cacc << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0; + + // Get the centroid for the feature + centroid[0] = centroids[3 * featureId + 0]; + centroid[1] = centroids[3 * featureId + 1]; + centroid[2] = centroids[3 * featureId + 2]; + + // for each triangle we need the transformation matrix A defined by the three points as columns + // Loop over all triangle faces + int32_t tCount = 0; + for(usize i = 0; i < numFaces; i++) + { + if(faceLabels[2 * i] != featureId && faceLabels[2 * i + 1] != featureId) + { + continue; + } + tCount++; + usize compIndex = (faceLabels[2 * i] == featureId ? 0 : 1); + std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); - const nx::core::Point3Df& a = vertCoords[0] - centroid; - const nx::core::Point3Df& b = (compIndex == 0 ? vertCoords[1] : vertCoords[2]) - centroid; - const nx::core::Point3Df& c = (compIndex == 0 ? vertCoords[2] : vertCoords[1]) - centroid; + const nx::core::Point3Df& a = vertCoords[0] - centroid; + const nx::core::Point3Df& b = (compIndex == 0 ? vertCoords[1] : vertCoords[2]) - centroid; + const nx::core::Point3Df& c = (compIndex == 0 ? vertCoords[2] : vertCoords[1]) - centroid; - Matrix3x3 A; - A << a[0], b[0], c[0], a[1], b[1], c[1], a[2], b[2], c[2]; + Matrix3x3 A; + A << a[0], b[0], c[0], a[1], b[1], c[1], a[2], b[2], c[2]; - float64 dA = A.determinant(); + float64 dA = A.determinant(); - Cacc = (Cacc + dA * (A * (C * (A.transpose())))).eval(); - Vol += (dA / 6.0f); // accumulate the volumes - } - - Cacc = (Cacc / Vol).eval(); - Matrix3x3 Cinertia = ID * Cacc.trace() - Cacc; - // extract the moments from the inertia tensor - Eigen::Vector3d e(Cinertia(0, 0), Cinertia(1, 1), Cinertia(2, 2)); - auto sols = CC * e; - omega3S[featureId] = static_cast(((Vol * Vol) / sols.prod()) / k_Sphere); + Cacc = (Cacc + dA * (A * (C * (A.transpose())))).eval(); + Vol += (dA / 6.0f); // accumulate the volumes + } - m_FeatureMoments[featureId * 6 + 0] = sols[0]; - m_FeatureMoments[featureId * 6 + 1] = sols[1]; - m_FeatureMoments[featureId * 6 + 2] = sols[2]; - m_FeatureMoments[featureId * 6 + 3] = -Cinertia(0, 1); - m_FeatureMoments[featureId * 6 + 4] = -Cinertia(0, 2); - m_FeatureMoments[featureId * 6 + 5] = -Cinertia(1, 2); + Cacc = (Cacc / Vol).eval(); + Cinertia = ID * Cacc.trace() - Cacc; + // extract the moments from the inertia tensor + Eigen::Vector3d e(Cinertia(0, 0), Cinertia(1, 1), Cinertia(2, 2)); + auto sols = CC * e; + omega3S[featureId] = static_cast(((Vol * Vol) / sols.prod()) / k_Sphere); + } /** * This next section finds the principle axis via eigenvalues. @@ -169,28 +241,22 @@ void ComputeTriangleGeomShapes::findMoments() * The main goal is to derive the eigenvalues from the moment of inertia tensor therein finding the eigenvectors, * which are the angular velocity vectors. * - * Code Keynotes for reviewers: - * - Hessenburg Decomposition is pre-processing to get an upper/right triangular matrix. - * - This significantly reduces the iterations needed for QR Decomposition - * - QR Factorization is different than QR Decomposition. - * - QR Decomposition expresses the product of an Orthogonal Matrix (Q) and an right/upper triangular matrix (R) as a singular matrix (A). + * Hessenburg Decomposition is pre-processing to get an upper/right triangular matrix. + * QR Decomposition expresses the product of an Orthogonal Matrix (Q) and an right/upper triangular matrix (R) as a singular matrix (A). */ + // !!! 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); - // TODO: - // - Remove Copying between processing (inline integration) - // - Extract real part of complexes stored in eigenvalues/vectors - Eigen::HessenbergDecomposition hessDecomp(Cinertia); - Matrix3x3 hessenMatrixUpper = hessDecomp.matrixH(); + // pre-processing to get an upper/right triangular matrix. + Eigen::HessenbergDecomposition hessDecomp(Cinertia); // Pass in Cinertia for implicit calculation // It is important to maintain the FullPivHouseholderQR implementation here, it is rank preserving and most numerically stable - Eigen::FullPivHouseholderQR hQR(hessenMatrixUpper); - Matrix3x3 hqrMatrix = hQR.matrixQR(); + Eigen::FullPivHouseholderQR hQR(hessDecomp.matrixH()); // pass in Hessenburg Matrix (Upper) for implicit calculation // Extract eigenvalues and eigenvectors - Eigen::EigenSolver eigenSolver(hqrMatrix); + Eigen::EigenSolver eigenSolver(hQR.matrixQR()); // pass in HQR Matrix for implicit calculation // The primary axis is the largest eigenvalue Eigen::EigenSolver::EigenvalueType eigenvalues = eigenSolver.eigenvalues(); @@ -227,243 +293,90 @@ void ComputeTriangleGeomShapes::findMoments() // double condNumNorms = hqrMatrix.completeOrthogonalDecomposition().pseudoInverse().norm() * hqrMatrix.norm(); // // std::cout << "\n Condition Number from Norms: " << condNumNorms << std::endl; - } -} -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findAxes() -{ - constexpr float64 k_ScaleFactor = 1.0; - const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); - - auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); - auto& aspectRatios = m_DataStructure.getDataRefAs(m_InputValues->AspectRatiosArrayPath); - - usize numFeatures = centroids.getNumberOfTuples(); - - m_FeatureEigenVals.resize(numFeatures * 3); + // 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); - for(usize i = 1; i < numFeatures; i++) - { - if(m_ShouldCancel) - { - return; - } - float64 ixx = m_FeatureMoments[i * 6 + 0]; - float64 iyy = m_FeatureMoments[i * 6 + 1]; - float64 izz = m_FeatureMoments[i * 6 + 2]; - - float64 ixy = m_FeatureMoments[i * 6 + 3]; - float64 iyz = m_FeatureMoments[i * 6 + 4]; - float64 ixz = m_FeatureMoments[i * 6 + 5]; - - float64 a = 1.0; - float64 b = (-ixx - iyy - izz); - float64 c = ((ixx * izz) + (ixx * iyy) + (iyy * izz) - (ixz * ixz) - (ixy * ixy) - (iyz * iyz)); - float64 d = ((ixz * iyy * ixz) + (ixy * izz * ixy) + (iyz * ixx * iyz) - (ixx * iyy * izz) - (ixy * iyz * ixz) - (ixy * iyz * ixz)); - - // f and g are the p and q values when reducing the cubic equation to t^3 + pt + q = 0 - float64 f = ((3.0 * c / a) - ((b / a) * (b / a))) / 3.0; - float64 g = ((2.0 * (b / a) * (b / a) * (b / a)) - (9.0 * b * c / (a * a)) + (27.0 * (d / a))) / 27.0; - float64 h = (g * g / 4.0) + (f * f * f / 27.0); - float64 rSquare = (g * g / 4.0) - h; - float64 r = sqrt(rSquare); - if(rSquare < 0.0) - { - r = 0.0; - } - float64 theta = 0; - if(r == 0) - { - theta = 0; - } - if(r != 0) + /** + * The following section finds axes + */ { - float64 value = -g / (2.0 * r); - if(value > 1) + if(m_ShouldCancel) { - value = 1.0; + return {}; } - if(value < -1) + // 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) { - value = -1.0; + bOverA = 0.0f, cOverA = 0.0f; } - theta = acos(value); - } - float64 const1 = pow(r, 0.33333333333); - float64 const2 = cos(theta / 3.0); - float64 const3 = b / (3.0 * a); - float64 const4 = 1.7320508 * sin(theta / 3.0); - - float64 r1 = 2 * const1 * const2 - (const3); - float64 r2 = -const1 * (const2 - (const4)) - const3; - float64 r3 = -const1 * (const2 + (const4)) - const3; - m_FeatureEigenVals[3 * i] = r1; - m_FeatureEigenVals[3 * i + 1] = r2; - m_FeatureEigenVals[3 * i + 2] = r3; - - float64 i1 = (15.0 * r1) / (4.0 * M_PI); - float64 i2 = (15.0 * r2) / (4.0 * M_PI); - float64 i3 = (15.0 * r3) / (4.0 * M_PI); - float64 aRatio = (i1 + i2 - i3) / 2.0f; - float64 bRatio = (i1 + i3 - i2) / 2.0f; - float64 cRatio = (i2 + i3 - i1) / 2.0f; - a = (aRatio * aRatio * aRatio * aRatio) / (bRatio * cRatio); - a = pow(a, 0.1); - b = bRatio / aRatio; - b = sqrt(b) * a; - c = aRatio / (a * a * a * b); - - axisLengths[3 * i] = static_cast(a / k_ScaleFactor); - axisLengths[3 * i + 1] = static_cast(b / k_ScaleFactor); - axisLengths[3 * i + 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; } - aspectRatios[2 * i] = bOverA; - aspectRatios[2 * i + 1] = cOverA; - } -} -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findAxisEulers() -{ - const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); - usize numFeatures = centroids.getNumberOfTuples(); - - auto& axisEulerAngles = m_DataStructure.getDataRefAs(m_InputValues->AxisEulerAnglesArrayPath); - - for(usize i = 1; i < numFeatures; i++) - { - if(m_ShouldCancel) - { - return; - } - float64 ixx = m_FeatureMoments[i * 6 + 0]; - float64 iyy = m_FeatureMoments[i * 6 + 1]; - float64 izz = m_FeatureMoments[i * 6 + 2]; - float64 ixy = m_FeatureMoments[i * 6 + 3]; - float64 iyz = m_FeatureMoments[i * 6 + 4]; - float64 ixz = m_FeatureMoments[i * 6 + 5]; - float64 radius1 = m_FeatureEigenVals[3 * i]; - float64 radius2 = m_FeatureEigenVals[3 * i + 1]; - float64 radius3 = m_FeatureEigenVals[3 * i + 2]; - - float64 e[3][1]; - float64 vect[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - e[0][0] = radius1; - e[1][0] = radius2; - e[2][0] = radius3; - float64 uber[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - float64 bMatrix[3]; - bMatrix[0] = 0.0000001; - bMatrix[1] = 0.0000001; - bMatrix[2] = 0.0000001; - - for(int32 j = 0; j < 3; j++) + /** + * The following section calculates the axis eulers + */ { - uber[0][0] = ixx - e[j][0]; - uber[0][1] = ixy; - uber[0][2] = ixz; - uber[1][0] = ixy; - uber[1][1] = iyy - e[j][0]; - uber[1][2] = iyz; - uber[2][0] = ixz; - uber[2][1] = iyz; - uber[2][2] = izz - e[j][0]; - std::array, 3> uberelim{}; - std::array, 3> uberbelim{}; - int32 elimCount = 0; - - for(int32 a = 0; a < 3; a++) + if(m_ShouldCancel) { - for(int32 b = 0; b < 3; b++) - { - uberelim[elimCount][b] = uber[a][b]; - } - uberbelim[elimCount][0] = bMatrix[a]; - elimCount++; - } - for(int32 k = 0; k < elimCount - 1; k++) - { - for(int32 l = k + 1; l < elimCount; l++) - { - float64 c = uberelim[l][k] / uberelim[k][k]; - for(int32 r = k + 1; r < elimCount; r++) - { - uberelim[l][r] = uberelim[l][r] - c * uberelim[k][r]; - } - uberbelim[l][0] = uberbelim[l][0] - c * uberbelim[k][0]; - } - } - uberbelim[elimCount - 1][0] = uberbelim[elimCount - 1][0] / uberelim[elimCount - 1][elimCount - 1]; - for(int32 l = 1; l < elimCount; l++) - { - int32 r = (elimCount - 1) - l; - float64 sum = 0.0; - for(int32 n = r + 1; n < elimCount; n++) - { - sum = sum + (uberelim[r][n] * uberbelim[n][0]); - } - uberbelim[r][0] = (uberbelim[r][0] - sum) / uberelim[r][r]; + return {}; } - for(int32 p = 0; p < elimCount; p++) + + // EigenVector associated with the largest EigenValue goes in the 3rd column + auto col3 = eigenvectors.col(idxs[0]); + + // Then the next largest into the 2nd column + auto col2 = eigenvectors.col(idxs[1]); + + // The smallest into the 1rst column + auto col1 = eigenvectors.col(idxs[2]); + + // 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()}}; + // clang-format on + + // check for right-handedness + OrientationTransformation::ResultType result = OrientationTransformation::om_check(OrientationD(g)); + if(result.result == 0) { - vect[j][p] = uberbelim[p][0]; + g[2][0] *= -1.0f; + g[2][1] *= -1.0f; + g[2][2] *= -1.0f; } - } - float64 n1X = vect[0][0]; - float64 n1Y = vect[0][1]; - float64 n1Z = vect[0][2]; - float64 n2X = vect[1][0]; - float64 n2Y = vect[1][1]; - float64 n2Z = vect[1][2]; - float64 n3X = vect[2][0]; - float64 n3Y = vect[2][1]; - float64 n3Z = vect[2][2]; - float64 norm1 = sqrt(((n1X * n1X) + (n1Y * n1Y) + (n1Z * n1Z))); - float64 norm2 = sqrt(((n2X * n2X) + (n2Y * n2Y) + (n2Z * n2Z))); - float64 norm3 = sqrt(((n3X * n3X) + (n3Y * n3Y) + (n3Z * n3Z))); - n1X = n1X / norm1; - n1Y = n1Y / norm1; - n1Z = n1Z / norm1; - n2X = n2X / norm2; - n2Y = n2Y / norm2; - n2Z = n2Z / norm2; - n3X = n3X / norm3; - n3Y = n3Y / norm3; - n3Z = n3Z / norm3; - - // 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) - float32 g[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - g[0][0] = static_cast(n3X); - g[0][1] = static_cast(n3Y); - g[0][2] = static_cast(n3Z); - g[1][0] = static_cast(n2X); - g[1][1] = static_cast(n2Y); - g[1][2] = static_cast(n2Z); - g[2][0] = static_cast(n1X); - g[2][1] = static_cast(n1Y); - g[2][2] = static_cast(n1Z); - - // check for right-handedness - OrientationTransformation::ResultType result = OrientationTransformation::om_check(OrientationF(g)); - if(result.result == 0) - { - g[2][0] *= -1.0f; - g[2][1] *= -1.0f; - g[2][2] *= -1.0f; - } - - auto euler = OrientationTransformation::om2eu(OrientationF(g)); + auto euler = OrientationTransformation::om2eu(OrientationD(g)); - axisEulerAngles[3 * i] = euler[0]; - axisEulerAngles[3 * i + 1] = euler[1]; - axisEulerAngles[3 * i + 2] = euler[2]; + axisEulerAngles[3 * featureId] = euler[0]; + axisEulerAngles[3 * featureId + 1] = euler[1]; + axisEulerAngles[3 * featureId + 2] = euler[2]; + } } + + return {}; } diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp index d48bd1f139..75cb133e1c 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp @@ -45,13 +45,5 @@ class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapes const ComputeTriangleGeomShapesInputValues* m_InputValues = nullptr; const std::atomic_bool& m_ShouldCancel; const IFilter::MessageHandler& m_MessageHandler; - - std::vector m_FeatureMoments; - std::vector m_FeatureEigenVals; - - void findMoments(); - void findAxes(); - void findAxisEulers(); }; - } // namespace nx::core From 9c4809241da92a3ff49b08028a53f027c91b5107 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Tue, 31 Dec 2024 17:03:59 -0500 Subject: [PATCH 07/12] strip unnecessary eigen processing --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 53 ++++++------------- .../ComputeTriangleGeomShapesFilter.cpp | 3 +- 2 files changed, 18 insertions(+), 38 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 65ae8e05f1..9afe29b2ed 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -240,23 +240,13 @@ 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. - * - * Hessenburg Decomposition is pre-processing to get an upper/right triangular matrix. - * QR Decomposition expresses the product of an Orthogonal Matrix (Q) and an right/upper triangular matrix (R) as a singular matrix (A). */ // !!! 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); - // pre-processing to get an upper/right triangular matrix. - Eigen::HessenbergDecomposition hessDecomp(Cinertia); // Pass in Cinertia for implicit calculation - - // It is important to maintain the FullPivHouseholderQR implementation here, it is rank preserving and most numerically stable - Eigen::FullPivHouseholderQR hQR(hessDecomp.matrixH()); // pass in Hessenburg Matrix (Upper) for implicit calculation - - // Extract eigenvalues and eigenvectors - Eigen::EigenSolver eigenSolver(hQR.matrixQR()); // pass in HQR Matrix for implicit calculation + Eigen::EigenSolver eigenSolver(Cinertia); // pass in HQR Matrix for implicit calculation // The primary axis is the largest eigenvalue Eigen::EigenSolver::EigenvalueType eigenvalues = eigenSolver.eigenvalues(); @@ -267,32 +257,21 @@ Result<> ComputeTriangleGeomShapes::operator()() /** * Following section for debugging */ - // 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; - // - // // Condition Number calculations from here: https://stackoverflow.com/questions/33575478/how-can-you-find-the-condition-number-in-eigen - // // Calculate Condition Number from JacobiSVD - // Eigen::JacobiSVD svd(hqrMatrix); - // double cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size() - 1); - // std::cout << "\n Condition Number from JacobiSVD: " << cond << std::endl; - // - // // Calculate Condition Number from Norms - // double condNumNorms = hqrMatrix.completeOrthogonalDecomposition().pseudoInverse().norm() * hqrMatrix.norm(); + // std::cout << "Eigenvalues:\n" << eigenvalues << std::endl; + // std::cout << "\n Eigenvectors:\n" << eigenvectors << std::endl; // - // std::cout << "\n Condition Number from Norms: " << condNumNorms << 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 @@ -306,7 +285,7 @@ Result<> ComputeTriangleGeomShapes::operator()() { return {}; } - // Formula: I = (15.0 * eigenvalue) / 4.0 * Pi; + // 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; diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp index 91cc7490b4..f061fe9999 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp @@ -5,6 +5,7 @@ #include "simplnx/DataStructure/DataPath.hpp" #include "simplnx/Filter/Actions/CreateArrayAction.hpp" #include "simplnx/Parameters/ArraySelectionParameter.hpp" +#include "simplnx/Parameters/AttributeMatrixSelectionParameter.hpp" #include "simplnx/Parameters/DataGroupSelectionParameter.hpp" #include "simplnx/Parameters/DataObjectNameParameter.hpp" #include "simplnx/Parameters/GeometrySelectionParameter.hpp" @@ -57,7 +58,7 @@ Parameters ComputeTriangleGeomShapesFilter::parameters() const params.insertSeparator(Parameters::Separator{"Input Face Feature Data"}); params.insert(std::make_unique( k_FeatureAttributeMatrixPath_Key, "Face Feature Attribute Matrix", "The DataPath to the AttributeMatrix that holds feature data for the faces", - DataPath({"TriangleDataContainer", "Face Feature Data"}), DataGroupSelectionParameter::AllowedTypes{BaseGroup::GroupType::AttributeMatrix})); + DataPath({"TriangleDataContainer", "Face Feature Data"}))); params.insert(std::make_unique(k_CentroidsArrayPath_Key, "Face Feature Centroids", "Input DataPath to the **Feature Centroids** for the face data", DataPath({"Face Feature Data", "Centroids"}), ArraySelectionParameter::AllowedTypes{DataType::float32})); params.insert(std::make_unique(k_VolumesArrayPath_Key, "Face Feature Volumes", "Input DataPath to the **Feature Volumes** for the face data", From 74792969698f7437424e23b14a49a68ccdc865e2 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Thu, 2 Jan 2025 11:00:25 -0500 Subject: [PATCH 08/12] Correct CAD shape finder, Axis lengths are actually un-normalized ratios --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 115 +++++++++--------- 1 file changed, 55 insertions(+), 60 deletions(-) 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 {}; From 3dee71a24a8fe1065ab893d5ff65383e253b6544 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Thu, 2 Jan 2025 11:12:31 -0500 Subject: [PATCH 09/12] parameter update: remove unused volumes array --- .../Algorithms/ComputeTriangleGeomShapes.hpp | 1 - .../ComputeTriangleGeomShapesFilter.cpp | 25 ++++++++----------- .../ComputeTriangleGeomShapesFilter.hpp | 1 - 3 files changed, 10 insertions(+), 17 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp index 75cb133e1c..21aecd0ca1 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp @@ -15,7 +15,6 @@ struct ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesInputValues DataPath FaceLabelsArrayPath; DataPath FeatureAttributeMatrixPath; DataPath CentroidsArrayPath; - DataPath VolumesArrayPath; DataPath Omega3sArrayPath; DataPath AxisLengthsArrayPath; DataPath AxisEulerAnglesArrayPath; diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp index f061fe9999..0742af4930 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp @@ -56,13 +56,11 @@ Parameters ComputeTriangleGeomShapesFilter::parameters() const ArraySelectionParameter::AllowedTypes{nx::core::DataType::int32}, ArraySelectionParameter::AllowedComponentShapes{{2}})); params.insertSeparator(Parameters::Separator{"Input Face Feature Data"}); - params.insert(std::make_unique( - k_FeatureAttributeMatrixPath_Key, "Face Feature Attribute Matrix", "The DataPath to the AttributeMatrix that holds feature data for the faces", - DataPath({"TriangleDataContainer", "Face Feature Data"}))); + params.insert(std::make_unique(k_FeatureAttributeMatrixPath_Key, "Face Feature Attribute Matrix", + "The DataPath to the AttributeMatrix that holds feature data for the faces", + DataPath({"TriangleDataContainer", "Face Feature Data"}))); params.insert(std::make_unique(k_CentroidsArrayPath_Key, "Face Feature Centroids", "Input DataPath to the **Feature Centroids** for the face data", DataPath({"Face Feature Data", "Centroids"}), ArraySelectionParameter::AllowedTypes{DataType::float32})); - params.insert(std::make_unique(k_VolumesArrayPath_Key, "Face Feature Volumes", "Input DataPath to the **Feature Volumes** for the face data", - DataPath({"Face Feature Data", "Volumes"}), ArraySelectionParameter::AllowedTypes{DataType::float32})); params.insertSeparator(Parameters::Separator{"Output Face Feature Data"}); params.insert(std::make_unique(k_Omega3sArrayName_Key, "Omega3s", "The name of the DataArray that holds the calculated Omega3 values", "Omega3s")); @@ -77,7 +75,11 @@ Parameters ComputeTriangleGeomShapesFilter::parameters() const //------------------------------------------------------------------------------ IFilter::VersionType ComputeTriangleGeomShapesFilter::parametersVersion() const { - return 1; + return 2; + + // Version 1 -> 2 + // Change 1: + // Removed input volumes array } //------------------------------------------------------------------------------ @@ -94,14 +96,11 @@ IFilter::PreflightResult ComputeTriangleGeomShapesFilter::preflightImpl(const Da auto pFaceLabelsArrayPathValue = filterArgs.value(k_FaceLabelsArrayPath_Key); auto pFeatureAttributeMatrixPath = filterArgs.value(k_FeatureAttributeMatrixPath_Key); auto pCentroidsArrayPathValue = filterArgs.value(k_CentroidsArrayPath_Key); - auto pVolumesArrayPathValue = filterArgs.value(k_VolumesArrayPath_Key); auto omega3sArrayNameValue = filterArgs.value(k_Omega3sArrayName_Key); auto axisLengthsArrayNameValue = filterArgs.value(k_AxisLengthsArrayName_Key); auto axisEulerAnglesArrayNameValue = filterArgs.value(k_AxisEulerAnglesArrayName_Key); auto aspectRatiosArrayNameValue = filterArgs.value(k_AspectRatiosArrayName_Key); - PreflightResult preflightResult; - nx::core::Result resultOutputActions; // Ensure the Face Feature Attribute Matrix is really an AttributeMatrix @@ -142,11 +141,8 @@ IFilter::PreflightResult ComputeTriangleGeomShapesFilter::preflightImpl(const Da resultOutputActions.value().appendAction(std::move(createArrayAction)); } - // No preflight updated values are generated in this filter - std::vector preflightUpdatedValues; - - // Return both the resultOutputActions and the preflightUpdatedValues via std::move() - return {std::move(resultOutputActions), std::move(preflightUpdatedValues)}; + // Return both the resultOutputActions via std::move() + return {std::move(resultOutputActions)}; } //------------------------------------------------------------------------------ @@ -158,7 +154,6 @@ Result<> ComputeTriangleGeomShapesFilter::executeImpl(DataStructure& dataStructu inputValues.FaceLabelsArrayPath = filterArgs.value(k_FaceLabelsArrayPath_Key); inputValues.FeatureAttributeMatrixPath = filterArgs.value(k_FeatureAttributeMatrixPath_Key); inputValues.CentroidsArrayPath = filterArgs.value(k_CentroidsArrayPath_Key); - inputValues.VolumesArrayPath = filterArgs.value(k_VolumesArrayPath_Key); auto omega3sArrayNameValue = filterArgs.value(k_Omega3sArrayName_Key); auto axisLengthsArrayNameValue = filterArgs.value(k_AxisLengthsArrayName_Key); diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp index fc9d090f4d..1d1bda319e 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp @@ -27,7 +27,6 @@ class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesFilter : public IFilte static inline constexpr StringLiteral k_FaceLabelsArrayPath_Key = "face_labels_array_path"; static inline constexpr StringLiteral k_FeatureAttributeMatrixPath_Key = "feature_attribute_matrix_path"; static inline constexpr StringLiteral k_CentroidsArrayPath_Key = "centroids_array_path"; - static inline constexpr StringLiteral k_VolumesArrayPath_Key = "volumes_array_path"; static inline constexpr StringLiteral k_Omega3sArrayName_Key = "omega3s_array_name"; static inline constexpr StringLiteral k_AxisLengthsArrayName_Key = "axis_lengths_array_name"; static inline constexpr StringLiteral k_AxisEulerAnglesArrayName_Key = "axis_euler_angles_array_name"; From e8a3457565c1ecb1b54d4b3feff95277ff6468fd Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Fri, 3 Jan 2025 16:50:01 -0500 Subject: [PATCH 10/12] Patch test and calculate Euler Characteristic fast watertight validation --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 145 +++++++++++++++++- .../ComputeTriangleGeomShapesFilter.cpp | 1 - .../test/ComputeTriangleGeomShapesTest.cpp | 3 +- 3 files changed, 144 insertions(+), 5 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 0e28f6ad3e..7dfb9f885a 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -19,6 +19,138 @@ using VertsStore = AbstractDataStore(std::numeric_limits::max()); + +template +usize SafeEdgeCount(const AbstractDataStore& faceStore) +{ + const usize numFaces = faceStore.getNumberOfTuples(); + const usize numComp = faceStore.getNumberOfComponents(); + T v0 = 0; + T v1 = 0; + + std::set> edgeSet; + + for(usize i = 0; i < numFaces; i++) + { + const usize offset = i * numComp; + + for(usize j = 0; j < numComp; j++) + { + if(j == (numComp - 1)) + { + if(faceStore[offset + j] > faceStore[offset + 0]) + { + v0 = faceStore[offset + 0]; + v1 = faceStore[offset + j]; + } + else + { + v0 = faceStore[offset + j]; + v1 = faceStore[offset + 0]; + } + } + else + { + if(faceStore[offset + j] > faceStore[offset + j + 1]) + { + v0 = faceStore[offset + j + 1]; + v1 = faceStore[offset + j]; + } + else + { + v0 = faceStore[offset + j]; + v1 = faceStore[offset + j + 1]; + } + } + std::pair edge = std::make_pair(v0, v1); + edgeSet.insert(edge); + } + } + + return edgeSet.size(); +} + +template +usize FastEdgeCount(const AbstractDataStore& faceStore) +{ + const usize numFaces = faceStore.getNumberOfTuples(); + const usize numComp = faceStore.getNumberOfComponents(); + uint32 v0 = 0; + uint32 v1 = 0; + + std::set edgeSet; + + for(usize i = 0; i < numFaces; i++) + { + const usize offset = i * numComp; + + for(usize j = 0; j < numComp; j++) + { + if(j == (numComp - 1)) + { + if(faceStore[offset + j] > faceStore[offset + 0]) + { + v0 = static_cast(faceStore[offset + 0]); + v1 = static_cast(faceStore[offset + j]); + } + else + { + v0 = static_cast(faceStore[offset + j]); + v1 = static_cast(faceStore[offset + 0]); + } + } + else + { + if(faceStore[offset + j] > faceStore[offset + j + 1]) + { + v0 = static_cast(faceStore[offset + j + 1]); + v1 = static_cast(faceStore[offset + j]); + } + else + { + v0 = static_cast(faceStore[offset + j]); + v1 = static_cast(faceStore[offset + j + 1]); + } + } + + edgeSet.insert(static_cast(v0) << 32 | v1); + } + } + + return edgeSet.size(); +} + +template +usize FindNumEdges(const AbstractDataStore& faceStore, usize numVertices = (k_MaxOptimizedValue + 1)) +{ + // This case may seem niche, but it is designed with Indexing types in mind specifically IGeometry::MeshIndexType + if constexpr(!std::is_signed_v) + { + if(numVertices < k_MaxOptimizedValue) + { + // speedier method because max vertices value fits into uint32 + return FastEdgeCount(faceStore); + } + } + // Slower method to avoid overflow + return SafeEdgeCount(faceStore); +} + +usize FindEulerCharacteristic(usize numVertices, usize numFaces, usize numRegions) +{ + return numVertices + numFaces - (2 * numRegions); +} + +template +bool ValidateMesh(const AbstractDataStore& faceStore, usize numVertices, usize numRegions) +{ + // Expensive call + usize numEdges = FindNumEdges(faceStore, numVertices); + + return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions); +} + /** * @brief This will extract the 3 vertices from a given triangle face of a triangle geometry. This is MUCH faster * than calling the function in the Triangle Geometry because of the dynamic_cast<> that goes on in that function. @@ -136,8 +268,17 @@ Result<> ComputeTriangleGeomShapes::operator()() const TriStore& triangleList = triangleGeom.getFacesRef().getDataStoreRef(); const VertsStore& verts = triangleGeom.getVerticesRef().getDataStoreRef(); - const auto& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath); - const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); + const auto& faceLabels = m_DataStructure.getDataRefAs(m_InputValues->FaceLabelsArrayPath).getDataStoreRef(); + const auto& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath).getDataStoreRef(); + + // the assumption here is face labels contains information on region ids, that it is contiguous in the values, and that 0 is an invalid id + // (ie the max function means that if the values in array are [1,2,4,5] it will assume there are 5 regions) + usize numRegions = *std::max_element(faceLabels.begin(), faceLabels.end()); + + if(!ValidateMesh(triangleList, verts.getNumberOfTuples(), numRegions)) + { + return MakeErrorResult(-64720, fmt::format("The Euler Characteristic of the shape was found to be unequal to 2, this implies the shape may not be watertight or is malformed.")); + } // Calculated Arrays auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp index 0742af4930..931a83cefe 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp @@ -92,7 +92,6 @@ IFilter::UniquePointer ComputeTriangleGeomShapesFilter::clone() const IFilter::PreflightResult ComputeTriangleGeomShapesFilter::preflightImpl(const DataStructure& dataStructure, const Arguments& filterArgs, const MessageHandler& messageHandler, const std::atomic_bool& shouldCancel) const { - auto pFaceLabelsArrayPathValue = filterArgs.value(k_FaceLabelsArrayPath_Key); auto pFeatureAttributeMatrixPath = filterArgs.value(k_FeatureAttributeMatrixPath_Key); auto pCentroidsArrayPathValue = filterArgs.value(k_CentroidsArrayPath_Key); diff --git a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp index 776e6659df..1de0dd199a 100644 --- a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp @@ -25,7 +25,6 @@ const std::string k_FaceLabelsName = "Face Labels"; const std::string k_FaceFeatureName = "Face Feature Data"; const std::string k_FaceDataName = "Face Data"; const std::string k_CentroidsArrayName = "Centroids"; -const std::string k_VolumesArrayName = "Volumes"; const std::string k_Omega3SArrayName = "Omega3s [NX Computed]"; const std::string k_AxisLengthsArrayName = "AxisLengths [NX Computed]"; @@ -88,7 +87,7 @@ TEST_CASE("OrientationAnalysis::ComputeTriangleGeomShapes", "[OrientationAnalysi args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_FeatureAttributeMatrixPath_Key, std::make_any(faceFeatureDataPath)); args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_CentroidsArrayPath_Key, std::make_any(faceFeatureDataPath.createChildPath(k_CentroidsArrayName))); - args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_VolumesArrayPath_Key, std::make_any(faceFeatureDataPath.createChildPath(k_VolumesArrayName))); + // Output Vars args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_Omega3sArrayName_Key, std::make_any(k_Omega3SArrayName)); args.insertOrAssign(ComputeTriangleGeomShapesFilter::k_AxisLengthsArrayName_Key, std::make_any(k_AxisLengthsArrayName)); From 33f62d487c3e180673ce91f723e19ed517f6790f Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Wed, 8 Jan 2025 16:10:41 -0500 Subject: [PATCH 11/12] use unordered set in fast edge counter, and move to utility file for shared use --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 127 +--------------- src/simplnx/Utilities/GeometryHelpers.hpp | 142 ++++++++++++++++++ 2 files changed, 148 insertions(+), 121 deletions(-) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 7dfb9f885a..0e5aef6c59 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -5,10 +5,13 @@ #include "simplnx/DataStructure/DataGroup.hpp" #include "simplnx/DataStructure/Geometry/IGeometry.hpp" #include "simplnx/DataStructure/Geometry/TriangleGeom.hpp" +#include "simplnx/Utilities/GeometryHelpers.hpp" #include "EbsdLib/Core/Orientation.hpp" #include "EbsdLib/Core/OrientationTransformation.hpp" +#include + using namespace nx::core; namespace @@ -19,134 +22,16 @@ using VertsStore = AbstractDataStore(std::numeric_limits::max()); - -template -usize SafeEdgeCount(const AbstractDataStore& faceStore) -{ - const usize numFaces = faceStore.getNumberOfTuples(); - const usize numComp = faceStore.getNumberOfComponents(); - T v0 = 0; - T v1 = 0; - - std::set> edgeSet; - - for(usize i = 0; i < numFaces; i++) - { - const usize offset = i * numComp; - - for(usize j = 0; j < numComp; j++) - { - if(j == (numComp - 1)) - { - if(faceStore[offset + j] > faceStore[offset + 0]) - { - v0 = faceStore[offset + 0]; - v1 = faceStore[offset + j]; - } - else - { - v0 = faceStore[offset + j]; - v1 = faceStore[offset + 0]; - } - } - else - { - if(faceStore[offset + j] > faceStore[offset + j + 1]) - { - v0 = faceStore[offset + j + 1]; - v1 = faceStore[offset + j]; - } - else - { - v0 = faceStore[offset + j]; - v1 = faceStore[offset + j + 1]; - } - } - std::pair edge = std::make_pair(v0, v1); - edgeSet.insert(edge); - } - } - - return edgeSet.size(); -} - -template -usize FastEdgeCount(const AbstractDataStore& faceStore) -{ - const usize numFaces = faceStore.getNumberOfTuples(); - const usize numComp = faceStore.getNumberOfComponents(); - uint32 v0 = 0; - uint32 v1 = 0; - - std::set edgeSet; - - for(usize i = 0; i < numFaces; i++) - { - const usize offset = i * numComp; - - for(usize j = 0; j < numComp; j++) - { - if(j == (numComp - 1)) - { - if(faceStore[offset + j] > faceStore[offset + 0]) - { - v0 = static_cast(faceStore[offset + 0]); - v1 = static_cast(faceStore[offset + j]); - } - else - { - v0 = static_cast(faceStore[offset + j]); - v1 = static_cast(faceStore[offset + 0]); - } - } - else - { - if(faceStore[offset + j] > faceStore[offset + j + 1]) - { - v0 = static_cast(faceStore[offset + j + 1]); - v1 = static_cast(faceStore[offset + j]); - } - else - { - v0 = static_cast(faceStore[offset + j]); - v1 = static_cast(faceStore[offset + j + 1]); - } - } - - edgeSet.insert(static_cast(v0) << 32 | v1); - } - } - - return edgeSet.size(); -} - -template -usize FindNumEdges(const AbstractDataStore& faceStore, usize numVertices = (k_MaxOptimizedValue + 1)) -{ - // This case may seem niche, but it is designed with Indexing types in mind specifically IGeometry::MeshIndexType - if constexpr(!std::is_signed_v) - { - if(numVertices < k_MaxOptimizedValue) - { - // speedier method because max vertices value fits into uint32 - return FastEdgeCount(faceStore); - } - } - // Slower method to avoid overflow - return SafeEdgeCount(faceStore); -} - usize FindEulerCharacteristic(usize numVertices, usize numFaces, usize numRegions) { return numVertices + numFaces - (2 * numRegions); } -template +template bool ValidateMesh(const AbstractDataStore& faceStore, usize numVertices, usize numRegions) { // Expensive call - usize numEdges = FindNumEdges(faceStore, numVertices); + usize numEdges = GeometryHelpers::Connectivity::FindNumEdges(faceStore, numVertices); return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions); } @@ -274,7 +159,7 @@ Result<> ComputeTriangleGeomShapes::operator()() // the assumption here is face labels contains information on region ids, that it is contiguous in the values, and that 0 is an invalid id // (ie the max function means that if the values in array are [1,2,4,5] it will assume there are 5 regions) usize numRegions = *std::max_element(faceLabels.begin(), faceLabels.end()); - + if(!ValidateMesh(triangleList, verts.getNumberOfTuples(), numRegions)) { return MakeErrorResult(-64720, fmt::format("The Euler Characteristic of the shape was found to be unequal to 2, this implies the shape may not be watertight or is malformed.")); diff --git a/src/simplnx/Utilities/GeometryHelpers.hpp b/src/simplnx/Utilities/GeometryHelpers.hpp index 0785087f4b..134a17e026 100644 --- a/src/simplnx/Utilities/GeometryHelpers.hpp +++ b/src/simplnx/Utilities/GeometryHelpers.hpp @@ -9,6 +9,8 @@ #include +#include + namespace nx::core { namespace GeometryHelpers @@ -32,6 +34,146 @@ SIMPLNX_EXPORT std::string GenerateGeometryInfo(const nx::core::SizeVec3& dims, namespace Connectivity { +namespace detail +{ +constexpr uint64 k_MaxOptimizedValue = static_cast(std::numeric_limits::max()); + +/** + * @brief !!! DO NOT CALL DIRECTLY !!! Prefer FindNumEdges() + * @tparam T indexing type of the face store + * @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList + * @return usize This is the number of edges + */ +template +usize SafeEdgeCount(const AbstractDataStore& faceStore) +{ + const usize numFaces = faceStore.getNumberOfTuples(); + const usize numComp = faceStore.getNumberOfComponents(); + T v0 = 0; + T v1 = 0; + + std::set> edgeSet; + + for(usize i = 0; i < numFaces; i++) + { + const usize offset = i * numComp; + + for(usize j = 0; j < numComp; j++) + { + if(j == (numComp - 1)) + { + if(faceStore[offset + j] > faceStore[offset + 0]) + { + v0 = faceStore[offset + 0]; + v1 = faceStore[offset + j]; + } + else + { + v0 = faceStore[offset + j]; + v1 = faceStore[offset + 0]; + } + } + else + { + if(faceStore[offset + j] > faceStore[offset + j + 1]) + { + v0 = faceStore[offset + j + 1]; + v1 = faceStore[offset + j]; + } + else + { + v0 = faceStore[offset + j]; + v1 = faceStore[offset + j + 1]; + } + } + std::pair edge = std::make_pair(v0, v1); + edgeSet.insert(edge); + } + } + + return edgeSet.size(); +} + +/** + * @brief !!! DO NOT CALL DIRECTLY !!! Prefer FindNumEdges() + * @tparam T indexing type of the face store + * @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList + * @return usize This is the number of edges + */ +template +usize FastEdgeCount(const AbstractDataStore& faceStore) +{ + const usize numFaces = faceStore.getNumberOfTuples(); + const usize numComp = faceStore.getNumberOfComponents(); + uint32 v0 = 0; + uint32 v1 = 0; + + std::unordered_set edgeSet; + + for(usize i = 0; i < numFaces; i++) + { + const usize offset = i * numComp; + + for(usize j = 0; j < numComp; j++) + { + if(j == (numComp - 1)) + { + if(faceStore[offset + j] > faceStore[offset + 0]) + { + v0 = static_cast(faceStore[offset + 0]); + v1 = static_cast(faceStore[offset + j]); + } + else + { + v0 = static_cast(faceStore[offset + j]); + v1 = static_cast(faceStore[offset + 0]); + } + } + else + { + if(faceStore[offset + j] > faceStore[offset + j + 1]) + { + v0 = static_cast(faceStore[offset + j + 1]); + v1 = static_cast(faceStore[offset + j]); + } + else + { + v0 = static_cast(faceStore[offset + j]); + v1 = static_cast(faceStore[offset + j + 1]); + } + } + + edgeSet.insert(static_cast(v0) << 32 | v1); + } + } + + return edgeSet.size(); +} +} // namespace detail + +/** + * @brief !!! EXPENSIVE !!! This function is a wrapper method for implicitly determining the correct edge calculation/counting algorithm + * @tparam T indexing type of the face store + * @param faceStore This is the face indexing list containing values that correspond to indices in the SharedVertexList + * @param numVertices This value is optional, since it is exclusively used for determining if faster algorithm is viable + * @return usize This is the number of edges + */ +template +usize FindNumEdges(const AbstractDataStore& faceStore, usize numVertices = (detail::k_MaxOptimizedValue + 1)) +{ + // This case may seem niche, but it is designed with Indexing types in mind specifically IGeometry::MeshIndexType + if constexpr(!std::is_signed_v) + { + if(numVertices < detail::k_MaxOptimizedValue) + { + // speedier method because max vertices value fits into uint32 + return detail::FastEdgeCount(faceStore); + } + } + // Slower method to avoid overflow + return detail::SafeEdgeCount(faceStore); +} + /** * @brief * @tparam T From f21b84d3d9a17ccc8a842ca248e4463fd03b2117 Mon Sep 17 00:00:00 2001 From: nyoungbq Date: Fri, 10 Jan 2025 17:15:20 -0500 Subject: [PATCH 12/12] [broken compile] intersection algorithm V1 --- .../Algorithms/ComputeTriangleGeomShapes.cpp | 160 ++++++++++++++++++ 1 file changed, 160 insertions(+) diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp index 0e5aef6c59..68507d0630 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -36,6 +36,166 @@ bool ValidateMesh(const AbstractDataStore& faceStore, usize numVertices, usiz return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions); } +// Be sure to pass validity bool +struct FeatureIntersectionTriangles +{ + IGeometry::MeshIndexType xIndex = 0; + IGeometry::MeshIndexType yIndex = 0; + IGeometry::MeshIndexType zIndex = 0; +}; + +std::array FindIntersections(const AbstractDataStore& faceLabelsStore, const AbstractDataStore& TriStore, + const AbstractDataStore& vertexStore, const AbstractDataStore& centroidsStore, + IGeometry::MeshIndexType featureId, std::vector& intersections) +{ + // Alias for templating later + using T = IGeometry::SharedVertexList::value_type; + + for(usize i = 0; i < faceLabelsStore.getNumberOfTuples(); i++) + { + // TODO: + // - Cancel Check Here + // - pass in orientation matrix + // - Sort out directional vectors + // - Potentially calc dist from feature centroid to triangle centroid + + if(faceLabelsStore[i] != featureId && faceLabelsStore[i + 1] != featureId) + { + // Triangle not in feature continue + continue; + } + + Eigen::Matrix orientationMatrix; + + // derive the direction vector for each corresponding axis + Eigen::Vector3d xDirVec = Eigen::Vector3d{1.0, 0.0, 0.0} * orientationMatrix; + Eigen::Vector3d yDirVec = Eigen::Vector3d{0.0, 1.0, 0.0} * orientationMatrix; + Eigen::Vector3d zDirVec = Eigen::Vector3d{0.0, 0.0, 1.0} * orientationMatrix; + + using PointT = Eigen::Vector3; + + PointT pointA = PointT{vertexStore[TriStore[i]], vertexStore[TriStore[i] + 1], vertexStore[TriStore[i] + 2]}; + PointT pointB = PointT{vertexStore[TriStore[i + 1]], vertexStore[TriStore[i + 1] + 1], vertexStore[TriStore[i + 1] + 2]}; + PointT pointC = PointT{vertexStore[TriStore[i + 2]], vertexStore[TriStore[i + 2] + 1], vertexStore[TriStore[i + 2] + 2]}; + + // Feature centroid + PointT origin = PointT{centroidsStore[i], centroidsStore[i + 1], centroidsStore[i + 2]}; + + constexpr T epsilon = std::numeric_limits::epsilon(); + + PointT edge1 = pointB - pointA; + PointT edge2 = pointC - pointA; + + PointT xCrossE2 = xDirVec.cross(edge2); + T xDet = edge1.dot(xCrossE2); + PointT yCrossE2 = yDirVec.cross(edge2); + T yDet = edge1.dot(yCrossE2); + PointT zCrossE2 = zDirVec.cross(edge2); + T zDet = edge1.dot(zCrossE2); + + if(xDet > -epsilon && xDet < epsilon) + { + // Ray is parallel to given triangle + return {}; + } + + if(yDet > -epsilon && yDet < epsilon) + { + // Ray is parallel to given triangle + return {}; + } + + if(zDet > -epsilon && zDet < epsilon) + { + // Ray is parallel to given triangle + return {}; + } + + PointT s = origin - pointA; + + T xInvDet = 1.0 / xDet; + T yInvDet = 1.0 / yDet; + T zInvDet = 1.0 / zDet; + + T xU = xInvDet * s.dot(xCrossE2); + T yU = yInvDet * s.dot(yCrossE2); + T zU = zInvDet * s.dot(zCrossE2); + + // Allow ADL for efficient absolute value function + using std::abs; + + if((xU < 0 && abs(xU) > epsilon) || (xU > 1 && abs(xU - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + if((yU < 0 && abs(yU) > epsilon) || (yU > 1 && abs(yU - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + if((zU < 0 && abs(zU) > epsilon) || (zU > 1 && abs(zU - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + PointT sCrossE1 = s.cross(edge1); + + T xV = xInvDet * xDirVec.dot(sCrossE1); + T yV = yInvDet * yDirVec.dot(sCrossE1); + T zV = zInvDet * zDirVec.dot(sCrossE1); + + if((xV < 0 && abs(xV) > epsilon) || (xU + xV > 1 && abs(xU + xV - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + if((yV < 0 && abs(yV) > epsilon) || (yU + yV > 1 && abs(yU + yV - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + if((zV < 0 && abs(zV) > epsilon) || (zU + zV > 1 && abs(zU + zV - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } + + T xT = xInvDet * edge2.dot(sCrossE1); + T yT = yInvDet * edge2.dot(sCrossE1); + T zT = zInvDet * edge2.dot(sCrossE1); + + if(xT > epsilon) + { + // Ray intersection + intersections[featureId].xIndex = i; + return origin + xDirVec * xT; + } + + if(yT > epsilon) + { + // Ray intersection + intersections[featureId].yIndex = i; + return origin + yDirVec * yT; + } + + if(zT > epsilon) + { + // Ray intersection + intersections[featureId].zIndex = i; + return origin + zDirVec * zT; + } + + // Line intersection but not ray intersection + return {}; + } +} + /** * @brief This will extract the 3 vertices from a given triangle face of a triangle geometry. This is MUCH faster * than calling the function in the Triangle Geometry because of the dynamic_cast<> that goes on in that function.