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..68507d0630 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp @@ -1,433 +1,291 @@ #include "ComputeTriangleGeomShapes.hpp" +#include "simplnx/Common/Constants.hpp" #include "simplnx/DataStructure/DataArray.hpp" #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 { +using TriStore = AbstractDataStore; +using VertsStore = AbstractDataStore; + +constexpr double k_Multiplier = 1.0 / (4.0 * Constants::k_PiD); 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; -// ----------------------------------------------------------------------------- +usize FindEulerCharacteristic(usize numVertices, usize numFaces, usize numRegions) +{ + return numVertices + numFaces - (2 * numRegions); +} + template -void FindTetrahedronInfo(const std::array& vertIds, const DataArray& vertPtr, const std::array& centroid, std::array& tetInfo) +bool ValidateMesh(const AbstractDataStore& faceStore, usize numVertices, usize numRegions) { - 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 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 + // Expensive call + usize numEdges = GeometryHelpers::Connectivity::FindNumEdges(faceStore, numVertices); - 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}; - - 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); - } + return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions); } -} // namespace - -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findMoments() +// Be sure to pass validity bool +struct FeatureIntersectionTriangles { - using MeshIndexType = IGeometry::MeshIndexType; - using SharedVertexListType = IGeometry::SharedVertexList; + 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; - 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); + 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 - auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); + if(faceLabelsStore[i] != featureId && faceLabelsStore[i + 1] != featureId) + { + // Triangle not in feature continue + continue; + } - usize numFaces = faceLabels.getNumberOfTuples(); - usize numFeatures = centroids.getNumberOfTuples(); - m_FeatureMoments.resize(numFeatures * 6, 0.0); + Eigen::Matrix orientationMatrix; - 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}; + // 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; - for(usize i = 0; i < numFaces; i++) - { - triangles.getFacePointIds(i, vertIds); - for(usize j = 0; j < 2; j++) + 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) { - if(j == 1) - { - std::swap(vertIds[2], vertIds[1]); - } - int32 gnum = faceLabels[2 * i + j]; - if(gnum > 0) - { - centroid[0] = centroids[3 * gnum + 0]; - centroid[1] = centroids[3 * gnum + 1]; - centroid[2] = centroids[3 * gnum + 2]; - FindTetrahedronInfo(vertIds, vertPtr, centroid, tetInfo); - for(usize iter = 0; iter < 8; iter++) - { - 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]); - - 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); - - 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]); - } - } + // Ray is parallel to given triangle + return {}; } - } - 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]; - 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 omega3 = vol5 / o3; - omega3 = omega3 / k_Sphere; - if(omega3 > 1) + if(yDet > -epsilon && yDet < epsilon) { - omega3 = 1.0; + // Ray is parallel to given triangle + return {}; } - if(vol5 == 0.0) + + if(zDet > -epsilon && zDet < epsilon) { - omega3 = 0.0; + // Ray is parallel to given triangle + return {}; } - omega3S[i] = static_cast(omega3); - } -} -// ----------------------------------------------------------------------------- -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; + PointT s = origin - pointA; - const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); + T xInvDet = 1.0 / xDet; + T yInvDet = 1.0 / yDet; + T zInvDet = 1.0 / zDet; - auto& axisLengths = m_DataStructure.getDataRefAs(m_InputValues->AxisLengthsArrayPath); - auto& aspectRatios = m_DataStructure.getDataRefAs(m_InputValues->AspectRatiosArrayPath); + T xU = xInvDet * s.dot(xCrossE2); + T yU = yInvDet * s.dot(yCrossE2); + T zU = zInvDet * s.dot(zCrossE2); - usize numFeatures = centroids.getNumberOfTuples(); + // Allow ADL for efficient absolute value function + using std::abs; - m_FeatureEigenVals.resize(numFeatures * 3); + if((xU < 0 && abs(xU) > epsilon) || (xU > 1 && abs(xU - 1.0) > epsilon)) + { + // Ray is parallel to given triangle + return {}; + } - for(usize i = 1; i < numFeatures; i++) - { - 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) + if((yU < 0 && abs(yU) > epsilon) || (yU > 1 && abs(yU - 1.0) > epsilon)) { - r = 0.0; + // Ray is parallel to given triangle + return {}; } - float64 theta = 0; - if(r == 0) + + if((zU < 0 && abs(zU) > epsilon) || (zU > 1 && abs(zU - 1.0) > epsilon)) { - theta = 0; + // Ray is parallel to given triangle + return {}; } - if(r != 0) + + 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)) { - float64 value = -g / (2.0 * r); - if(value > 1) - { - value = 1.0; - } - if(value < -1) - { - value = -1.0; - } - theta = acos(value); + // 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; } - 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) + + if(yT > epsilon) { - bOverA = 0.0f, cOverA = 0.0f; + // Ray intersection + intersections[featureId].yIndex = i; + return origin + yDirVec * yT; } - aspectRatios[2 * i] = bOverA; - aspectRatios[2 * i + 1] = cOverA; + + if(zT > epsilon) + { + // Ray intersection + intersections[featureId].zIndex = i; + return origin + zDirVec * zT; + } + + // Line intersection but not ray intersection + return {}; } } -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findAxisEulers() +/** + * @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) { - const Float32Array& centroids = m_DataStructure.getDataRefAs(m_InputValues->CentroidsArrayPath); - usize numFeatures = centroids.getNumberOfTuples(); - - auto& axisEulerAngles = m_DataStructure.getDataRefAs(m_InputValues->AxisEulerAnglesArrayPath); + 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]}}; +} - for(usize i = 1; i < numFeatures; i++) +/** + * @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) { - 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++) + // sorted[2] = a; + if(b > c) { - 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++) - { - 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]; - } - for(int32 p = 0; p < elimCount; p++) - { - vect[j][p] = uberbelim[p][0]; - } + // sorted[1] = b; + // sorted[0] = c; + idx = {C, B, A}; } - - 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) + else { - g[2][0] *= -1.0f; - g[2][1] *= -1.0f; - g[2][2] *= -1.0f; + // 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}; + } - auto euler = OrientationTransformation::om2eu(OrientationF(g)); - - axisEulerAngles[3 * i] = euler[0]; - axisEulerAngles[3 * i + 1] = euler[1]; - axisEulerAngles[3 * i + 2] = euler[2]; + if(!lowToHigh) + { + std::swap(idx[0], idx[2]); } + return idx; } +} // namespace + // ----------------------------------------------------------------------------- ComputeTriangleGeomShapes::ComputeTriangleGeomShapes(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ComputeTriangleGeomShapesInputValues* inputValues) @@ -450,9 +308,235 @@ const std::atomic_bool& ComputeTriangleGeomShapes::getCancel() // ----------------------------------------------------------------------------- Result<> ComputeTriangleGeomShapes::operator()() { - findMoments(); - findAxes(); - findAxisEulers(); + using MeshIndexType = IGeometry::MeshIndexType; + 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).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); + 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(); + + nx::core::Point3Df centroid = {0.0F, 0.0F, 0.0F}; + + // 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; + + // 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 "Feature" which is the number of tuples in the "Centroids" array + // We could parallelize over the features? + for(usize featureId = 1; featureId < numFeatures; featureId++) + { + /** + * The following section calculates moment of inertia tensor (Cinertia) and omega3s + */ + { + 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++) + { + 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; + + 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(); + + Cacc = (Cacc + dA * (A * (C * (A.transpose())))).eval(); + Vol += (dA / 6.0f); // accumulate the volumes + } + + 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. + * 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. + */ + Eigen::EigenSolver eigenSolver(Cinertia); + + // 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(); + + /** + * 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; + + // 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 calculates the axis eulers + */ + { + if(m_ShouldCancel) + { + return {}; + } + + // 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) + { + g[2][0] *= -1.0f; + g[2][1] *= -1.0f; + g[2][2] *= -1.0f; + } + + auto euler = OrientationTransformation::om2eu(OrientationD(g)); + + axisEulerAngles[3 * featureId] = euler[0]; + 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 {}; } diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp index 839e83cb53..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; @@ -45,12 +44,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 diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp index 91cc7490b4..931a83cefe 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" @@ -55,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"}), DataGroupSelectionParameter::AllowedTypes{BaseGroup::GroupType::AttributeMatrix})); + 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")); @@ -76,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 } //------------------------------------------------------------------------------ @@ -89,18 +92,14 @@ 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); - 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 @@ -141,11 +140,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)}; } //------------------------------------------------------------------------------ @@ -157,7 +153,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"; 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..1de0dd199a 100644 --- a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp @@ -8,28 +8,54 @@ #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"; const std::string k_Omega3SArrayName = "Omega3s [NX Computed]"; 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 +63,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))); + + // 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 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 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