diff --git a/CMakeLists.txt b/CMakeLists.txt index 3f827ff4d5..bef254c41e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -529,6 +529,7 @@ set(SIMPLNX_HDRS ${SIMPLNX_SOURCE_DIR}/Utilities/HistogramUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/MemoryUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/StringUtilities.hpp + ${SIMPLNX_SOURCE_DIR}/Utilities/IntersectionUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/IParallelAlgorithm.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/ParallelDataAlgorithm.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/ParallelData2DAlgorithm.hpp @@ -548,7 +549,6 @@ set(SIMPLNX_HDRS ${SIMPLNX_SOURCE_DIR}/Utilities/ClusteringUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/MontageUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/SIMPLConversion.hpp - ${SIMPLNX_SOURCE_DIR}/Utilities/IntersectionUtilities.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/Math/GeometryMath.hpp ${SIMPLNX_SOURCE_DIR}/Utilities/Math/MatrixMath.hpp diff --git a/src/Plugins/OrientationAnalysis/CMakeLists.txt b/src/Plugins/OrientationAnalysis/CMakeLists.txt index 52962201ad..54df4c050b 100644 --- a/src/Plugins/OrientationAnalysis/CMakeLists.txt +++ b/src/Plugins/OrientationAnalysis/CMakeLists.txt @@ -52,7 +52,7 @@ set(FilterList ComputeSchmidsFilter ComputeShapesFilter ComputeSlipTransmissionMetricsFilter - # ComputeTriangleGeomShapesFilter + ComputeShapesTriangleGeomFilter ConvertHexGridToSquareGridFilter ConvertOrientationsFilter ConvertQuaternionFilter @@ -186,7 +186,7 @@ set(filter_algorithms ComputeSchmids ComputeShapes ComputeSlipTransmissionMetrics - # ComputeTriangleGeomShapes + ComputeShapesTriangleGeom ConvertHexGridToSquareGrid ConvertQuaternion CreateEnsembleInfo diff --git a/src/Plugins/OrientationAnalysis/docs/ComputeShapesTriangleGeomFilter.md b/src/Plugins/OrientationAnalysis/docs/ComputeShapesTriangleGeomFilter.md new file mode 100644 index 0000000000..c3ca9977ee --- /dev/null +++ b/src/Plugins/OrientationAnalysis/docs/ComputeShapesTriangleGeomFilter.md @@ -0,0 +1,88 @@ +# Compute Feature Shapes from Triangle Geometry + +## Group (Subgroup) + +Statistics (Morphological) + +## Warning + +This filter has two caveats. + +Firstly, the axial lengths of this filter will be different than those produced by the voxelized version of this filter. This is for two reasons: + +- The sampling rate and density for the grid that was used to voxelize the mesh. See *Sample Triangle Geometry on Regular Grid* (RegularGridSampleSurfaceMesh). +- This filter determines axial lengths via distance from feature centroid to mesh intersection points along each of the principle axes. This means they are relative to the mesh itself rather than the grid it exists in. + +Secondly, shapes that exhibit rotational symmetry (e.g. cube, sphere, regular octahedron, etc.) may have different Euler Angles than those of the voxelized implementation, but they are functionally identical. This is more prevalent in meshes with less traingles, but this is seemly due to the fact the tested shapes are more uniform in low-poly. It is presumed that fiducial markers will stabilize ouputs for these specific shapes. + +## Description + +This **Filter** calculates the second-order moments of each enclosed **Feature** in a **Triangle Geometry**. The +second-order moments allow for the determination of the *principal axis lengths, principal axis directions, aspect +ratios and moment invariant Omega3s*. The *principal axis lengths* are those of a "best-fit" ellipsoid. The algorithm +for determining the moments and these values is as follows: + +1. For each **Triangle** on the bounding surface of a **Feature**, construct a tetrahedron whose fourth vertex is the + centroid of the **Feature**, ensuring normals are consistent (this **Filter** uses the convention where normals point + inwards; note that the actual winding of the **Triangle Geometry** is not modified) +2. Subdivide each constructed tetrahedron into 8 smaller tetrahedra +3. For each subdivided tetrahedron, compute the distance from that tetrahedron's centroid to the centroid of the + parent **Feature** +4. For each subdivided tetrahedron, calculate Ixx, Iyy, Izz, Ixy, Ixz and Iyz using the x, y and z distances determined + in step 1 +5. Use the relationship of *principal moments* to the *principal axis lengths* for an ellipsoid, which can be found + in [4], to determine the *Axis Lengths* +6. Calculate the *Aspect Ratios* from the *Axis Lengths* found in step 5. +7. Determine the Euler angles required to represent the *principal axis directions* in the *sample reference frame* and + store them as the **Feature**'s *Axis Euler Angles*. +8. Calculate the moment invariant Omega3 as defined in [2] and is discussed further in [1] and [3] + +*Note:* Due to the method used to subdivide the tetrahedra, some sharp corners of shapes may not be properly +represented, resulting in inaccurate Omega3 values. This problem is especially apparent for perfect rectangular prisms, +but any shape with clear sharp corners may be affected. + +% Auto generated parameter table will be inserted here + +## References + +[1] Representation and Reconstruction of Three-dimensional Microstructures in Ni-based Superalloys, AFOSR +FA9550-07-1-0179 Final Report, 20 Dec 2010. + +[2] On the use of moment invariants for the automated classification of 3-D particle shapes, J. MacSleyne, J.P. Simmons +and M. De Graef, Modeling and Simulations in Materials Science and Engineering, 16, 045008 (2008). + +[3] n-Dimensional Moment Invariants and Conceptual Mathematical Theory of Recognition n-Dimensional Solids, Alexander G. +Mamistvalov, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 20, NO. 8, AUGUST 1998, p. 819-831. + +[4] M. Groeber, M. Uchic, D. Dimiduk, and S. Ghosh. A Framework for Automated Analysis and Simulation of 3D +Polycrystalline Microstructures, Part 1: Statistical Characterization Acta Materialia, 56 (2008), 1257-1273. + +## Example Pipelines + +## Exemplars For Test/Validation + +Note that the `AxisEulerAngles` are in radians. The `_NUMBER` is the component index. The names correspond to files in the test data. The test data also contains a `.dream3d` file (and corresponding `.xdmf` file) that has the preprocessed input data as well as two vertex geometries that correspond to the axial intersection points for both this filter (named `stl`) and the voxelized implementation (named `voxelized`) to demonstrate the difference between the two results. See `validation` folder in the associated test file. The code to generate this was not production level so it was not included in the release of this filter. + +| Index | Exemplar Omega3s | Exemplar AxisEulerAngles_0 | Exemplar AxisEulerAngles_1 | Exemplar AxisEulerAngles_2 | Exemplar AxisLengths_0 | Exemplar AxisLengths_1 | Exemplar AxisLengths_2 | Centroids_0 | Centroids_1 | Centroids_2 | STL File List | +|---|---|---|---|---|---|---|---|---|---|---|--------| +| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 'JUNK' | +| 1 | 0.78787351 | 6.2519455 | 0.46384287 | 1.5707964 | 0.55903065 | 0.5 | 0.5590716 | 11 | 1.5 | 0.5 | '101_cube.stl' | +| 2 | 0.90795749 | 1.5712752 | 1.5711994 | 4.711997 | 0.50000107 | 0.50000107 | 0.50000536 | 7.9999948 | 1.500001 | 0.50000101 | '102_rounded_cube.stl' | +| 3 | 0.99998772 | 3.1415451 | 1.5707964 | 1.5707963 | 0.4999963 | 0.50000113 | 0.49999821 | 5.0000029 | 1.4999999 | 0.50000352 | '103_sphere.stl' | +| 4 | 0.81056947 | 0 | 0 | 0 | 0.5 | 0.5 | 0.5 | 2 | 1.5 | 0.5 | '104_octahedron.stl' | +| 5 | 0.78787351 | 0 | 1.5707964 | 4.712389 | 1.5 | 1 | 0.5 | 11 | 3.5 | 1.5 | '105_rectangular_prism.stl' | +| 6 | 0.90794492 | 2.3219979e-10 | 1.5707964 | 1.5707964 | 1.5000002 | 1.0000105 | 0.50000501 | 8.0000105 | 3.500005 | 1.4999998 | '106_rounded_cube_elongated.stl' | +| 7 | 0.99998951 | 3.1415927 | 1.5707964 | 1.5707964 | 1.499998 | 0.99999982 | 0.50000066 | 5 | 3.4999993 | 1.5000019 | '107_ellipsoid.stl' | +| 8 | 0.4052847 | 2.9784336 | 1.2776501 | 2.7590237 | 0.72980016 | 0.63307023 | 0.90308428 | 11.949251 | 5.9832501 | 0.56225002 | '108_tetrahedron.stl' | +| 9 | 0.86399871 | 4.7129421 | 1.5707964 | 1.5707964 | 1.25 | 1.0001171 | 1.0003718 | 2.0074062 | 6.1888757 | 1.25 | '109_cylinder.stl' | +| 10 | 0.96753073 | 4.712389 | 1.7409153 | 4.712389 | 1.0468618 | 1.1038886 | 1.2814574 | 8.2501669 | 6.1889997 | 1.03175 | '110_icosahedron.stl' | +| 11 | 0.95603704 | 1.5717114 | 2.5809844 | 4.7131243 | 0.95104223 | 0.95071417 | 0.94695413 | 4.8930006 | 6.1890001 | 0.80900002 | '111_dodecahedron.stl' | +| 12 | 0.81056947 | 3.1415927 | 1.5707964 | 1.5707964 | 1.5 | 1 | 0.5 | 2 | 3.5 | 1.5 | '112_octahedron.stl' | + +## License & Copyright + +Please see the description file distributed with this plugin. + +## DREAM3D-NX Help + +If you need help, need to file a bug report or want to request a new feature, please head over to the [DREAM3DNX-Issues](https://github.com/BlueQuartzSoftware/DREAM3DNX-Issues/discussions) GitHub site where the community of DREAM3D-NX users can help answer your questions. diff --git a/src/Plugins/OrientationAnalysis/docs/ComputeTriangleGeomShapesFilter.md b/src/Plugins/OrientationAnalysis/docs/ComputeTriangleGeomShapesFilter.md deleted file mode 100644 index 8e00622fb2..0000000000 --- a/src/Plugins/OrientationAnalysis/docs/ComputeTriangleGeomShapesFilter.md +++ /dev/null @@ -1,57 +0,0 @@ -# Compute Feature Shapes from Triangle Geometry - -## Group (Subgroup) - -Statistics (Morphological) - -## Description - -This **Filter** calculates the second-order moments of each enclosed **Feature** in a **Triangle Geometry**. The -second-order moments allow for the determination of the *principal axis lengths, principal axis directions, aspect -ratios and moment invariant Omega3s*. The *principal axis lengths* are those of a "best-fit" ellipsoid. The algorithm -for determining the moments and these values is as follows: - -1. For each **Triangle** on the bounding surface of a **Feature**, construct a tetrahedron whose fourth vertex is the - centroid of the **Feature**, ensuring normals are consistent (this **Filter** uses the convention where normals point - inwards; note that the actual winding of the **Triangle Geometry** is not modified) -2. Subdivide each constructed tetrahedron into 8 smaller tetrahedra -3. For each subdivided tetrahedron, compute the distance from that tetrahedron's centroid to the centroid of the - parent **Feature** -4. For each subdivided tetrahedron, calculate Ixx, Iyy, Izz, Ixy, Ixz and Iyz using the x, y and z distances determined - in step 1 -5. Use the relationship of *principal moments* to the *principal axis lengths* for an ellipsoid, which can be found - in [4], to determine the *Axis Lengths* -6. Calculate the *Aspect Ratios* from the *Axis Lengths* found in step 5. -7. Determine the Euler angles required to represent the *principal axis directions* in the *sample reference frame* and - store them as the **Feature**'s *Axis Euler Angles*. -8. Calculate the moment invariant Omega3 as defined in [2] and is discussed further in [1] and [3] - -*Note:* Due to the method used to subdivide the tetrahedra, some sharp corners of shapes may not be properly -represented, resulting in inaccurate Omega3 values. This problem is especially apparent for perfect rectangular prisms, -but any shape with clear sharp corners may be affected. - -% Auto generated parameter table will be inserted here - -## References ## - -[1] Representation and Reconstruction of Three-dimensional Microstructures in Ni-based Superalloys, AFOSR -FA9550-07-1-0179 Final Report, 20 Dec 2010. - -[2] On the use of moment invariants for the automated classification of 3-D particle shapes, J. MacSleyne, J.P. Simmons -and M. De Graef, Modeling and Simulations in Materials Science and Engineering, 16, 045008 (2008). - -[3] n-Dimensional Moment Invariants and Conceptual Mathematical Theory of Recognition n-Dimensional Solids, Alexander G. -Mamistvalov, IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. 20, NO. 8, AUGUST 1998, p. 819-831. - -[4] M. Groeber, M. Uchic, D. Dimiduk, and S. Ghosh. A Framework for Automated Analysis and Simulation of 3D -Polycrystalline Microstructures, Part 1: Statistical Characterization Acta Materialia, 56 (2008), 1257-1273. - -## Example Pipelines - -## License & Copyright - -Please see the description file distributed with this plugin. - -## DREAM3D-NX Help - -If you need help, need to file a bug report or want to request a new feature, please head over to the [DREAM3DNX-Issues](https://github.com/BlueQuartzSoftware/DREAM3DNX-Issues/discussions) GitHub site where the community of DREAM3D-NX users can help answer your questions. 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/ComputeShapesTriangleGeom.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapesTriangleGeom.cpp new file mode 100644 index 0000000000..6a14ec736e --- /dev/null +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapesTriangleGeom.cpp @@ -0,0 +1,458 @@ +#include "ComputeShapesTriangleGeom.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 "simplnx/Utilities/IntersectionUtilities.hpp" + +#include "EbsdLib/Core/Orientation.hpp" +#include "EbsdLib/Core/OrientationTransformation.hpp" + +#include + +using namespace nx::core; + +namespace +{ +using TriStore = AbstractDataStore; +using VertsStore = AbstractDataStore; + +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 + const usize numEdges = GeometryHelpers::Connectivity::FindNumEdges(faceStore, numVertices); + + return numEdges == FindEulerCharacteristic(numVertices, faceStore.getNumberOfTuples(), numRegions); +} + +struct AxialLengths +{ + IGeometry::SharedVertexList::value_type xLength = 0.0; + IGeometry::SharedVertexList::value_type yLength = 0.0; + IGeometry::SharedVertexList::value_type zLength = 0.0; +}; + +// Eigen implementation of Moller-Trumbore intersection algorithm adapted to account for distance +template +AxialLengths FindIntersections(const Eigen::Matrix& orientationMatrix, const AbstractDataStore& faceLabelsStore, + const AbstractDataStore& triStore, const AbstractDataStore& vertexStore, + const AbstractDataStore& centroidsStore, IGeometry::MeshIndexType featureId, const std::atomic_bool& shouldCancel) +{ + constexpr T epsilon = std::numeric_limits::epsilon(); + + AxialLengths lengths; + + using PointT = Eigen::Vector3; + + // derive the direction vector for each corresponding axis (using unit vectors) + const PointT xDirVec = PointT{1.0, 0.0, 0.0}.transpose() * orientationMatrix; + const PointT yDirVec = PointT{0.0, 1.0, 0.0}.transpose() * orientationMatrix; + const PointT zDirVec = PointT{0.0, 0.0, 1.0}.transpose() * orientationMatrix; + + IntersectionUtilities::MTPointsCache cache; + + // Feature Centroid + cache.origin = PointT{centroidsStore[3 * featureId], centroidsStore[(3 * featureId) + 1], centroidsStore[(3 * featureId) + 2]}; + + for(usize i = 0; i < faceLabelsStore.getNumberOfTuples(); i++) + { + if(shouldCancel) + { + return lengths; + } + + if(faceLabelsStore[2 * i] != featureId && faceLabelsStore[(2 * i) + 1] != featureId) + { + // Triangle not in feature continue + continue; + } + + // Here we are manually extracting the vertex points from the SharedVertexList + const usize threeCompIndex = 3 * i; + + // Extract tuple index of the first vertex and compute the index to the first x-value in the SharedVertexList + const usize vertAIndex = triStore[threeCompIndex] * 3; + cache.pointA = PointT{vertexStore[vertAIndex], vertexStore[vertAIndex + 1], vertexStore[vertAIndex + 2]}; + + const usize vertBIndex = triStore[threeCompIndex + 1] * 3; + cache.pointB = PointT{vertexStore[vertBIndex], vertexStore[vertBIndex + 1], vertexStore[vertBIndex + 2]}; + + const usize vertCIndex = triStore[threeCompIndex + 2] * 3; + cache.pointC = PointT{vertexStore[vertCIndex], vertexStore[vertCIndex + 1], vertexStore[vertCIndex + 2]}; + + cache.edge1 = cache.pointB - cache.pointA; + cache.edge2 = cache.pointC - cache.pointA; + + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + const std::optional xIntersect = IntersectionUtilities::MTIntersection(xDirVec, cache); + if(xIntersect.has_value()) + { + // Ray intersection + T distance = std::sqrt((xIntersect.value() - cache.origin).array().square().sum()); + if(abs(distance) > std::abs(lengths.xLength)) + { + lengths.xLength = distance; + } + } + + const std::optional yIntersect = IntersectionUtilities::MTIntersection(yDirVec, cache); + if(yIntersect.has_value()) + { + // Ray intersection + T distance = std::sqrt((yIntersect.value() - cache.origin).array().square().sum()); + if(abs(distance) > std::abs(lengths.yLength)) + { + lengths.yLength = distance; + } + } + + const std::optional zIntersect = IntersectionUtilities::MTIntersection(zDirVec, cache); + if(zIntersect.has_value()) + { + // Ray intersection + T distance = std::sqrt((zIntersect.value() - cache.origin).array().square().sum()); + if(abs(distance) > std::abs(lengths.zLength)) + { + lengths.zLength = distance; + } + } + } + + return lengths; +} + +/** + * @brief This will extract the 3 vertices from aVal given triangle face of aVal 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 usize v0Idx = triangleList[triangleId * 3]; + const usize v1Idx = triangleList[(triangleId * 3) + 1]; + const 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]}}; +} + +/** + * @brief Sorts the 3 values + * @param aVal First Value + * @param bVal Second Value + * @param cVal Third Value + * @return The indices in their sorted order + */ +template +std::array TripletSort(T aVal, T bVal, T cVal, bool lowToHigh) +{ + constexpr size_t k_AIdx = 0; + constexpr size_t k_BIdx = 1; + constexpr size_t k_CIdx = 2; + std::array idx = {0, 1, 2}; + if(aVal > bVal && aVal > cVal) + { + // sorted[2] = aVal; + if(bVal > cVal) + { + // sorted[1] = bVal; + // sorted[0] = cVal; + idx = {k_CIdx, k_BIdx, k_AIdx}; + } + else + { + // sorted[1] = cVal; + // sorted[0] = bVal; + idx = {k_BIdx, k_CIdx, k_AIdx}; + } + } + else if(bVal > aVal && bVal > cVal) + { + // sorted[2] = bVal; + if(aVal > cVal) + { + // sorted[1] = aVal; + // sorted[0] = cVal; + idx = {k_CIdx, k_AIdx, k_BIdx}; + } + else + { + // sorted[1] = cVal; + // sorted[0] = aVal; + idx = {k_AIdx, k_CIdx, k_BIdx}; + } + } + else if(aVal > bVal) + { + // sorted[1] = aVal; + // sorted[0] = bVal; + // sorted[2] = cVal; + idx = {k_BIdx, k_AIdx, k_CIdx}; + } + else if(aVal >= cVal && bVal >= cVal) + { + // sorted[0] = cVal; + // sorted[1] = aVal; + // sorted[2] = bVal; + idx = {k_CIdx, k_AIdx, k_BIdx}; + } + else + { + // sorted[0] = aVal; + // sorted[1] = bVal; + // sorted[2] = cVal; + idx = {k_AIdx, k_BIdx, k_CIdx}; + } + + if(!lowToHigh) + { + std::swap(idx[0], idx[2]); + } + return idx; +} + +} // namespace + +// ----------------------------------------------------------------------------- +ComputeShapesTriangleGeom::ComputeShapesTriangleGeom(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, + ComputeShapesTriangleGeomInputValues* inputValues) +: m_DataStructure(dataStructure) +, m_InputValues(inputValues) +, m_ShouldCancel(shouldCancel) +, m_MessageHandler(mesgHandler) +{ +} + +// ----------------------------------------------------------------------------- +ComputeShapesTriangleGeom::~ComputeShapesTriangleGeom() noexcept = default; + +// ----------------------------------------------------------------------------- +const std::atomic_bool& ComputeShapesTriangleGeom::getCancel() +{ + return m_ShouldCancel; +} + +// ----------------------------------------------------------------------------- +Result<> ComputeShapesTriangleGeom::operator()() +{ + 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) + const 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; + + const usize numFaces = faceLabels.getNumberOfTuples(); + const 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; + constexpr float64 k_Sphere = (2000.0 * M_PI * M_PI) / 9.0; + + // define the canonical cMatrix matrix + constexpr float64 aVal = 1.0 / 60.0; + constexpr float64 bVal = aVal / 2.0; + // clang-format off + Matrix3x3 cMatrix; + cMatrix << aVal, bVal, bVal, bVal, aVal, bVal, bVal, bVal, aVal; + + // and the identity matrix + Matrix3x3 identityMat; + identityMat << 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0; + + // The cMatrix-Prime matrix + Matrix3x3 cPrime; + cPrime << -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 {}; + } + float64 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++; + const usize compIndex = (faceLabels[2 * i] == featureId ? 0 : 1); + std::array vertCoords = GetFaceCoordinates(i, verts, triangleList); + + const nx::core::Point3Df& aVert = vertCoords[0] - centroid; + const nx::core::Point3Df& bVert = (compIndex == 0 ? vertCoords[1] : vertCoords[2]) - centroid; + const nx::core::Point3Df& cVert = (compIndex == 0 ? vertCoords[2] : vertCoords[1]) - centroid; + + Matrix3x3 aMat; + aMat << aVert[0], bVert[0], cVert[0], aVert[1], bVert[1], cVert[1], aVert[2], bVert[2], cVert[2]; + + const float64 detA = aMat.determinant(); + + Cacc = (Cacc + detA * (aMat * (cMatrix * (aMat.transpose())))).eval(); + Vol += (detA / 6.0f); // accumulate the volumes + } + + Cacc = (Cacc / Vol).eval(); + Cinertia = identityMat * Cacc.trace() - Cacc; + // extract the moments from the inertia tensor + const Eigen::Vector3d eVec(Cinertia(0, 0), Cinertia(1, 1), Cinertia(2, 2)); + auto sols = cPrime * eVec; + 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. + */ + const 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'; + // float64 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); + + Matrix3x3 orientationMatrix = {}; + + /** + * 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) + orientationMatrix.row(0) = col1.real(); + orientationMatrix.row(1) = col2.real(); + orientationMatrix.row(2) = col3.real(); + + auto euler = OrientationTransformation::om2eu(OrientationD(orientationMatrix.data(), 9)); + + axisEulerAngles[3 * featureId] = static_cast(euler[0]); + axisEulerAngles[(3 * featureId) + 1] = static_cast(euler[1]); + axisEulerAngles[(3 * featureId) + 2] = static_cast(euler[2]); + } + + /** + * The following section finds axes + */ + { + if(m_ShouldCancel) + { + return {}; + } + + const ::AxialLengths lengths = FindIntersections(orientationMatrix, faceLabels, triangleList, verts, centroids, featureId, m_ShouldCancel); + + // Check for zeroes (zeroes = probably invalid) + if(lengths.xLength == 0.0 || lengths.yLength == 0.0 || lengths.zLength == 0.0) + { + return MakeErrorResult(-64721, fmt::format("{}({}): One or more of the axis lengths for feature {} was unable to be found. This indicates the geometry was malformed.\nFeature Centroid(XYZ): " + "[{},{},{}]\nX Length: {}\nY Length: {}\nZ Length: {}", + __FILE__, __LINE__, featureId, centroids[(3 * featureId) + 0], centroids[(3 * featureId) + 1], centroids[(3 * featureId) + 2], lengths.xLength, + lengths.yLength, lengths.zLength)); + } + + axisLengths[3 * featureId] = static_cast(lengths.xLength); + axisLengths[(3 * featureId) + 1] = static_cast(lengths.yLength); + axisLengths[(3 * featureId) + 2] = static_cast(lengths.zLength); + auto bOverA = static_cast(lengths.yLength / lengths.xLength); + auto cOverA = static_cast(lengths.zLength / lengths.xLength); + 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/ComputeShapesTriangleGeom.hpp similarity index 52% rename from src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp rename to src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapesTriangleGeom.hpp index 839e83cb53..a6a56c557c 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeShapesTriangleGeom.hpp @@ -9,13 +9,12 @@ namespace nx::core { -struct ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesInputValues +struct ORIENTATIONANALYSIS_EXPORT ComputeShapesTriangleGeomInputValues { DataPath TriangleGeometryPath; DataPath FaceLabelsArrayPath; DataPath FeatureAttributeMatrixPath; DataPath CentroidsArrayPath; - DataPath VolumesArrayPath; DataPath Omega3sArrayPath; DataPath AxisLengthsArrayPath; DataPath AxisEulerAnglesArrayPath; @@ -25,16 +24,16 @@ struct ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesInputValues /** * @class */ -class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapes +class ORIENTATIONANALYSIS_EXPORT ComputeShapesTriangleGeom { public: - ComputeTriangleGeomShapes(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ComputeTriangleGeomShapesInputValues* inputValues); - ~ComputeTriangleGeomShapes() noexcept; + ComputeShapesTriangleGeom(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ComputeShapesTriangleGeomInputValues* inputValues); + ~ComputeShapesTriangleGeom() noexcept; - ComputeTriangleGeomShapes(const ComputeTriangleGeomShapes&) = delete; - ComputeTriangleGeomShapes(ComputeTriangleGeomShapes&&) noexcept = delete; - ComputeTriangleGeomShapes& operator=(const ComputeTriangleGeomShapes&) = delete; - ComputeTriangleGeomShapes& operator=(ComputeTriangleGeomShapes&&) noexcept = delete; + ComputeShapesTriangleGeom(const ComputeShapesTriangleGeom&) = delete; + ComputeShapesTriangleGeom(ComputeShapesTriangleGeom&&) noexcept = delete; + ComputeShapesTriangleGeom& operator=(const ComputeShapesTriangleGeom&) = delete; + ComputeShapesTriangleGeom& operator=(ComputeShapesTriangleGeom&&) noexcept = delete; Result<> operator()(); @@ -42,15 +41,8 @@ class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapes private: DataStructure& m_DataStructure; - const ComputeTriangleGeomShapesInputValues* m_InputValues = nullptr; + const ComputeShapesTriangleGeomInputValues* 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/Algorithms/ComputeTriangleGeomShapes.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp deleted file mode 100644 index 9211f382b5..0000000000 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.cpp +++ /dev/null @@ -1,458 +0,0 @@ -#include "ComputeTriangleGeomShapes.hpp" - -#include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/DataGroup.hpp" -#include "simplnx/DataStructure/Geometry/IGeometry.hpp" -#include "simplnx/DataStructure/Geometry/TriangleGeom.hpp" - -#include "EbsdLib/Core/Orientation.hpp" -#include "EbsdLib/Core/OrientationTransformation.hpp" - -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; - -// ----------------------------------------------------------------------------- -template -void FindTetrahedronInfo(const std::array& vertIds, const DataArray& vertPtr, 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 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 - - 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); - } -} - -} // namespace - -// ----------------------------------------------------------------------------- -void ComputeTriangleGeomShapes::findMoments() -{ - using MeshIndexType = IGeometry::MeshIndexType; - 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); - - auto& omega3S = m_DataStructure.getDataRefAs(m_InputValues->Omega3sArrayPath); - - usize numFaces = faceLabels.getNumberOfTuples(); - 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}; - - for(usize i = 0; i < numFaces; i++) - { - triangles.getFacePointIds(i, vertIds); - for(usize j = 0; j < 2; j++) - { - 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]); - } - } - } - } - - 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) - { - omega3 = 1.0; - } - if(vol5 == 0.0) - { - omega3 = 0.0; - } - 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; - - 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); - - 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) - { - r = 0.0; - } - float64 theta = 0; - if(r == 0) - { - theta = 0; - } - if(r != 0) - { - float64 value = -g / (2.0 * r); - if(value > 1) - { - value = 1.0; - } - if(value < -1) - { - value = -1.0; - } - 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 * 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++) - { - 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++) - { - 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]; - } - } - - 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)); - - axisEulerAngles[3 * i] = euler[0]; - axisEulerAngles[3 * i + 1] = euler[1]; - axisEulerAngles[3 * i + 2] = euler[2]; - } -} - -// ----------------------------------------------------------------------------- -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/ComputeShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesFilter.cpp index 2680e9eaef..7877e6db94 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesFilter.cpp @@ -37,13 +37,13 @@ Uuid ComputeShapesFilter::uuid() const //------------------------------------------------------------------------------ std::string ComputeShapesFilter::humanName() const { - return "Compute Feature Shapes"; + return "Compute Feature Shapes (Image Geometry)"; } //------------------------------------------------------------------------------ std::vector ComputeShapesFilter::defaultTags() const { - return {className(), "Statistics", "Morphological", "Find", "Generate", "Calculate", "Determine"}; + return {className(), "Statistics", "Morphological", "Find", "Generate", "Calculate", "Determine", "Omega3", "Axis Length"}; } //------------------------------------------------------------------------------ @@ -52,9 +52,9 @@ Parameters ComputeShapesFilter::parameters() const Parameters params; // Create the parameter descriptors that are needed for this filter - params.insertSeparator(Parameters::Separator{"Input Cell Data"}); params.insert(std::make_unique(k_SelectedImageGeometryPath_Key, "Selected Image Geometry", "The target geometry", DataPath{}, GeometrySelectionParameter::AllowedTypes{IGeometry::Type::Image})); + params.insertSeparator(Parameters::Separator{"Input Cell Data"}); params.insert(std::make_unique(k_CellFeatureIdsArrayPath_Key, "Cell Feature Ids", "Specifies to which Feature each Cell belongs", DataPath({"FeatureIds"}), ArraySelectionParameter::AllowedTypes{DataType::int32})); diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.cpp similarity index 78% rename from src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp rename to src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.cpp index 91cc7490b4..e9fc0e5b17 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.cpp @@ -1,10 +1,11 @@ -#include "ComputeTriangleGeomShapesFilter.hpp" -#include "OrientationAnalysis/Filters/Algorithms/ComputeTriangleGeomShapes.hpp" +#include "ComputeShapesTriangleGeomFilter.hpp" +#include "OrientationAnalysis/Filters/Algorithms/ComputeShapesTriangleGeom.hpp" #include "simplnx/DataStructure/AttributeMatrix.hpp" #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" @@ -14,37 +15,37 @@ using namespace nx::core; namespace nx::core { //------------------------------------------------------------------------------ -std::string ComputeTriangleGeomShapesFilter::name() const +std::string ComputeShapesTriangleGeomFilter::name() const { - return FilterTraits::name.str(); + return FilterTraits::name.str(); } //------------------------------------------------------------------------------ -std::string ComputeTriangleGeomShapesFilter::className() const +std::string ComputeShapesTriangleGeomFilter::className() const { - return FilterTraits::className; + return FilterTraits::className; } //------------------------------------------------------------------------------ -Uuid ComputeTriangleGeomShapesFilter::uuid() const +Uuid ComputeShapesTriangleGeomFilter::uuid() const { - return FilterTraits::uuid; + return FilterTraits::uuid; } //------------------------------------------------------------------------------ -std::string ComputeTriangleGeomShapesFilter::humanName() const +std::string ComputeShapesTriangleGeomFilter::humanName() const { - return "Compute Feature Shapes from Triangle Geometry"; + return "Compute Feature Shapes (Triangle Geometry)"; } //------------------------------------------------------------------------------ -std::vector ComputeTriangleGeomShapesFilter::defaultTags() const +std::vector ComputeShapesTriangleGeomFilter::defaultTags() const { - return {className(), "Statistics", "Morphological", "SurfaceMesh", "Find"}; + return {className(), "Statistics", "Morphological", "Find", "Generate", "Calculate", "Determine", "Omega3", "Axis Length", "Surface Mesh"}; } //------------------------------------------------------------------------------ -Parameters ComputeTriangleGeomShapesFilter::parameters() const +Parameters ComputeShapesTriangleGeomFilter::parameters() const { Parameters params; // Create the parameter descriptors that are needed for this filter @@ -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")); @@ -74,33 +73,33 @@ Parameters ComputeTriangleGeomShapesFilter::parameters() const } //------------------------------------------------------------------------------ -IFilter::VersionType ComputeTriangleGeomShapesFilter::parametersVersion() const +IFilter::VersionType ComputeShapesTriangleGeomFilter::parametersVersion() const { - return 1; + return 2; + + // Version 1 -> 2 + // Change 1: + // Removed input volumes array } //------------------------------------------------------------------------------ -IFilter::UniquePointer ComputeTriangleGeomShapesFilter::clone() const +IFilter::UniquePointer ComputeShapesTriangleGeomFilter::clone() const { - return std::make_unique(); + return std::make_unique(); } //------------------------------------------------------------------------------ -IFilter::PreflightResult ComputeTriangleGeomShapesFilter::preflightImpl(const DataStructure& dataStructure, const Arguments& filterArgs, const MessageHandler& messageHandler, +IFilter::PreflightResult ComputeShapesTriangleGeomFilter::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,23 +140,19 @@ 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)}; } //------------------------------------------------------------------------------ -Result<> ComputeTriangleGeomShapesFilter::executeImpl(DataStructure& dataStructure, const Arguments& filterArgs, const PipelineFilter* pipelineNode, const MessageHandler& messageHandler, +Result<> ComputeShapesTriangleGeomFilter::executeImpl(DataStructure& dataStructure, const Arguments& filterArgs, const PipelineFilter* pipelineNode, const MessageHandler& messageHandler, const std::atomic_bool& shouldCancel) const { - ComputeTriangleGeomShapesInputValues inputValues; + ComputeShapesTriangleGeomInputValues inputValues; inputValues.TriangleGeometryPath = filterArgs.value(k_TriGeometryDataPath_Key); 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); @@ -169,6 +164,6 @@ Result<> ComputeTriangleGeomShapesFilter::executeImpl(DataStructure& dataStructu inputValues.AxisEulerAnglesArrayPath = inputValues.FeatureAttributeMatrixPath.createChildPath(axisEulerAnglesArrayNameValue); inputValues.AspectRatiosArrayPath = inputValues.FeatureAttributeMatrixPath.createChildPath(aspectRatiosArrayNameValue); - return ComputeTriangleGeomShapes(dataStructure, messageHandler, shouldCancel, &inputValues)(); + return ComputeShapesTriangleGeom(dataStructure, messageHandler, shouldCancel, &inputValues)(); } } // namespace nx::core diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.hpp similarity index 84% rename from src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp rename to src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.hpp index fc9d090f4d..f854d5ee94 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.hpp @@ -8,26 +8,25 @@ namespace nx::core { /** - * @class ComputeTriangleGeomShapesFilter + * @class ComputeShapesTriangleGeomFilter * @brief This filter will .... */ -class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesFilter : public IFilter +class ORIENTATIONANALYSIS_EXPORT ComputeShapesTriangleGeomFilter : public IFilter { public: - ComputeTriangleGeomShapesFilter() = default; - ~ComputeTriangleGeomShapesFilter() noexcept override = default; + ComputeShapesTriangleGeomFilter() = default; + ~ComputeShapesTriangleGeomFilter() noexcept override = default; - ComputeTriangleGeomShapesFilter(const ComputeTriangleGeomShapesFilter&) = delete; - ComputeTriangleGeomShapesFilter(ComputeTriangleGeomShapesFilter&&) noexcept = delete; + ComputeShapesTriangleGeomFilter(const ComputeShapesTriangleGeomFilter&) = delete; + ComputeShapesTriangleGeomFilter(ComputeShapesTriangleGeomFilter&&) noexcept = delete; - ComputeTriangleGeomShapesFilter& operator=(const ComputeTriangleGeomShapesFilter&) = delete; - ComputeTriangleGeomShapesFilter& operator=(ComputeTriangleGeomShapesFilter&&) noexcept = delete; + ComputeShapesTriangleGeomFilter& operator=(const ComputeShapesTriangleGeomFilter&) = delete; + ComputeShapesTriangleGeomFilter& operator=(ComputeShapesTriangleGeomFilter&&) noexcept = delete; // Parameter Keys 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"; @@ -109,4 +108,4 @@ class ORIENTATIONANALYSIS_EXPORT ComputeTriangleGeomShapesFilter : public IFilte }; } // namespace nx::core -SIMPLNX_DEF_FILTER_TRAITS(nx::core, ComputeTriangleGeomShapesFilter, "e8f0fed3-d0d8-456e-b5a1-7961cc17b739"); +SIMPLNX_DEF_FILTER_TRAITS(nx::core, ComputeShapesTriangleGeomFilter, "e8f0fed3-d0d8-456e-b5a1-7961cc17b739"); diff --git a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt index aa13c0a38d..c3f7e12e63 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 + ComputeShapesTriangleGeomTest.cpp ConvertHexGridToSquareGridTest.cpp ConvertOrientationsTest.cpp ConvertQuaternionTest.cpp @@ -144,6 +144,7 @@ if(EXISTS "${DREAM3D_DATA_DIR}" AND SIMPLNX_DOWNLOAD_TEST_FILES) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME write_stats_gen_odf_angle_file.tar.gz SHA512 be3f663aae1f78e5b789200421534ed9fe293187ec3514796ac8177128b34ded18bb9a98b8e838bb283f9818ac30dc4b19ec379bdd581b1a98eb36d967cdd319) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME GrainMapper3D_Test_Files.tar.gz SHA512 cbd9cde2528d62de6cffecd125795b571f8b1a41c8eebbbbac5908f5382ac3c3c3ea56efc34b6e6ca713055f4fd13485aa6dd5acb9a69e9d4b3b16941ed1f96f) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_ReadChannel5_Test.tar.gz SHA512 c43533e8fc4c8c2f5a649a15fda856a023db1cbe972b82c34042ba050a25e7cacffd7cca979504f72a2384c7d7fd66308132f94366e367129d675f35abc82cba) + download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_compute_triangle_shapes_test.tar.gz SHA512 562ec65d22771817655ac85e267b1ea9822ab478ef732a52ce05052e4a171fc131bb879c9982df32b726d68a350f092695af97633f69e97226a6147cd0608b82) endif() diff --git a/src/Plugins/OrientationAnalysis/test/ComputeShapesTriangleGeomTest.cpp b/src/Plugins/OrientationAnalysis/test/ComputeShapesTriangleGeomTest.cpp new file mode 100644 index 0000000000..bf2765d98f --- /dev/null +++ b/src/Plugins/OrientationAnalysis/test/ComputeShapesTriangleGeomTest.cpp @@ -0,0 +1,95 @@ +#include "OrientationAnalysis/Filters/ComputeShapesTriangleGeomFilter.hpp" +#include "OrientationAnalysis/OrientationAnalysis_test_dirs.hpp" + +#include "simplnx/Parameters/ArrayCreationParameter.hpp" +#include "simplnx/Parameters/DataObjectNameParameter.hpp" +#include "simplnx/Parameters/GeometrySelectionParameter.hpp" +#include "simplnx/UnitTest/UnitTestCommon.hpp" + +#include + +#include +#include +#include + +namespace fs = std::filesystem; + +using namespace nx::core; + +#define SIMPLNX_WRITE_TEST_OUTPUT + +namespace ComputeShapesTriangleGeomFilterTest +{ +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_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 std::string k_ExemplarOmega3SArrayName = "Exemplar Omega3s"; +const std::string k_ExemplarAxisLengthsArrayName = "Exemplar AxisLengths"; +const std::string k_ExemplarAxisEulerAnglesArrayName = "Exemplar AxisEulerAngles"; +const std::string k_ExemplarAspectRatiosArrayName = "Exemplar AspectRatios"; + +constexpr StringLiteral k_GeomName = "InputGeometry"; + +const DataPath k_GeometryPath({k_GeomName}); + +const DataPath k_FaceFeatureAttributeMatrixPath = k_GeometryPath.createChildPath(k_FaceFeatureName); +const DataPath k_FaceDataPath = k_GeometryPath.createChildPath(k_FaceDataName); +const DataPath k_FaceLabelsPath = k_FaceDataPath.createChildPath(k_FaceLabelsName); +const DataPath k_FaceFeatureCentroidsPath = k_FaceFeatureAttributeMatrixPath.createChildPath(k_CentroidsArrayName); +} // namespace ComputeShapesTriangleGeomFilterTest + +using namespace ComputeShapesTriangleGeomFilterTest; + +// !!! See filter documentation for information on included data and how it was generated and visually validated !!! +TEST_CASE("OrientationAnalysis::ComputeShapesTriangleGeom", "[OrientationAnalysis][ComputeShapesTriangleGeom]") +{ + const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_CMakeExecutable, nx::core::unit_test::k_TestFilesDir, "7_compute_triangle_shapes_test.tar.gz", + "7_compute_triangle_shapes_test"); + + DataStructure exemplarDataStructure = UnitTest::LoadDataStructure(fs::path(fmt::format("{}/7_compute_triangle_shapes_test/test/7_exemplar_triangle_shapes.dream3d", unit_test::k_TestFilesDir))); + + // Instantiate the filter and an Arguments Object + ComputeShapesTriangleGeomFilter filter; + Arguments args; + + // Create default Parameters for the filter. + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_TriGeometryDataPath_Key, std::make_any(k_GeometryPath)); + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_FaceLabelsArrayPath_Key, std::make_any(k_FaceLabelsPath)); + + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_FeatureAttributeMatrixPath_Key, std::make_any(k_FaceFeatureAttributeMatrixPath)); + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_CentroidsArrayPath_Key, std::make_any(k_FaceFeatureCentroidsPath)); + + // Output Vars + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_Omega3sArrayName_Key, std::make_any(k_Omega3SArrayName)); + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_AxisLengthsArrayName_Key, std::make_any(k_AxisLengthsArrayName)); + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_AxisEulerAnglesArrayName_Key, std::make_any(k_AxisEulerAnglesArrayName)); + args.insertOrAssign(ComputeShapesTriangleGeomFilter::k_AspectRatiosArrayName_Key, std::make_any(k_AspectRatiosArrayName)); + + // Preflight the filter and check result + auto preflightResult = filter.preflight(exemplarDataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); + + // Execute the filter and check the result + auto executeResult = filter.execute(exemplarDataStructure, args); + SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result); + +#ifdef SIMPLNX_WRITE_TEST_OUTPUT + UnitTest::WriteTestDataStructure(exemplarDataStructure, fs::path(fmt::format("{}/{}.dream3d", unit_test::k_BinaryTestOutputDir, "ComputeShapesTriangleGeomTestOutput"))); +#endif + + UnitTest::CompareArrays(exemplarDataStructure, k_FaceFeatureAttributeMatrixPath.createChildPath(k_Omega3SArrayName), + k_FaceFeatureAttributeMatrixPath.createChildPath(k_ExemplarOmega3SArrayName)); + UnitTest::CompareArrays(exemplarDataStructure, k_FaceFeatureAttributeMatrixPath.createChildPath(k_AxisLengthsArrayName), + k_FaceFeatureAttributeMatrixPath.createChildPath(k_ExemplarAxisLengthsArrayName)); + UnitTest::CompareArrays(exemplarDataStructure, k_FaceFeatureAttributeMatrixPath.createChildPath(k_AxisEulerAnglesArrayName), + k_FaceFeatureAttributeMatrixPath.createChildPath(k_ExemplarAxisEulerAnglesArrayName)); + UnitTest::CompareArrays(exemplarDataStructure, k_FaceFeatureAttributeMatrixPath.createChildPath(k_AspectRatiosArrayName), + k_FaceFeatureAttributeMatrixPath.createChildPath(k_ExemplarAspectRatiosArrayName)); +} diff --git a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp b/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp deleted file mode 100644 index 473efe0c44..0000000000 --- a/src/Plugins/OrientationAnalysis/test/ComputeTriangleGeomShapesTest.cpp +++ /dev/null @@ -1,90 +0,0 @@ -#include "OrientationAnalysis/Filters/ComputeTriangleGeomShapesFilter.hpp" -#include "OrientationAnalysis/OrientationAnalysis_test_dirs.hpp" - -#include "simplnx/Parameters/ArrayCreationParameter.hpp" -#include "simplnx/Parameters/DataObjectNameParameter.hpp" -#include "simplnx/Parameters/GeometrySelectionParameter.hpp" -#include "simplnx/UnitTest/UnitTestCommon.hpp" - -#include - -using namespace nx::core; -using namespace nx::core::UnitTest; - -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_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); -} // namespace ComputeTriangleGeomShapesFilterTest - -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); - - { - // 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); - } - -#ifdef SIMPLNX_WRITE_TEST_OUTPUT - WriteTestDataStructure(dataStructure, fs::path(fmt::format("{}/find_triangle_geom_shapes.dream3d", unit_test::k_BinaryTestOutputDir))); -#endif -} 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..63e2b808f6 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 +{ +inline 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 d3adc0041d..1a6ed0728e 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 diff --git a/src/simplnx/Utilities/IntersectionUtilities.hpp b/src/simplnx/Utilities/IntersectionUtilities.hpp index 661723a083..ee91e94bfe 100644 --- a/src/simplnx/Utilities/IntersectionUtilities.hpp +++ b/src/simplnx/Utilities/IntersectionUtilities.hpp @@ -4,7 +4,10 @@ #include "simplnx/DataStructure/Geometry/IGeometry.hpp" #include "simplnx/DataStructure/Geometry/TriangleGeom.hpp" +#include + #include + namespace nx::core { namespace IntersectionUtilities @@ -217,6 +220,136 @@ inline bool RayTriangleIntersect2(const Vec3f& orig, const Vec3f& dir, const Vec return true; } +// These are values that are not associated with the direction of the ray +// Precalculating is only useful when calculating multiple rays (same origin/different directions) on the same triangle +template +struct MTPointsCache +{ + using PointT = Eigen::Vector3; + + // Ray origin + PointT origin; + + // Triangle Vertices + PointT pointA; + PointT pointB; + PointT pointC; + + // Specific edges + // edge1 = pointB - pointA; + PointT edge1; + // edge2 = pointC - pointA; + PointT edge2; + + // Pre-calculations + // sDist = origin - pointA; + PointT sDist; + // sCrossE1 = sDist.cross(edge1); + PointT sCrossE1; +}; + +// Eigen implementation of Moller-Trumbore intersection algorithm adapted to account for distance (Cached) +template +std::optional> MTIntersection(const Eigen::Vector3& dirVec, const MTPointsCache& cache) +{ + using PointT = Eigen::Vector3; + constexpr T epsilon = std::numeric_limits::epsilon(); + + PointT crossE2 = dirVec.cross(cache.edge2); + T determinant = cache.edge1.dot(crossE2); + + if(determinant > -epsilon && determinant < epsilon) + { + // Ray is parallel to given triangle + return {}; + } + + T invDet = 1.0 / determinant; + T uCoff = invDet * cache.sDist.dot(crossE2); + + // Allow ADL for efficient absolute value function + using std::abs; + if((uCoff < 0 && abs(uCoff) > epsilon) || (uCoff > 1 && abs(uCoff - 1.0) > epsilon)) + { + // Ray intersects plane, but not triangle + return {}; + } + + T vCoff = invDet * dirVec.dot(cache.sCrossE1); + + if((vCoff < 0 && abs(vCoff) > epsilon) || (uCoff + vCoff > 1 && abs(uCoff + vCoff - 1.0) > epsilon)) + { + // Ray intersects plane, but not triangle + return {}; + } + + T tCoff = invDet * cache.edge2.dot(cache.sCrossE1); + + if(tCoff > epsilon) + { + // Ray intersection + return {cache.origin + dirVec * tCoff}; + } + + // line intersection not ray intersection + return {}; +} + +// Eigen implementation of Moller-Trumbore intersection algorithm adapted to account for distance +template +std::optional> MTIntersection(const Eigen::Vector3& origin, const Eigen::Vector3& dirVec, const Eigen::Vector3& pointA, const Eigen::Vector3& pointB, + const Eigen::Vector3& pointC) +{ + using PointT = Eigen::Vector3; + constexpr T epsilon = std::numeric_limits::epsilon(); + + PointT edge1 = pointB - pointA; + PointT edge2 = pointC - pointA; + + PointT crossE2 = dirVec.cross(edge2); + T determinant = edge1.dot(crossE2); + + if(determinant > -epsilon && determinant < epsilon) + { + // Ray is parallel to given triangle + return {}; + } + + PointT sDist = origin - pointA; + + T invDet = 1.0 / determinant; + T uCoff = invDet * sDist.dot(crossE2); + + // Allow ADL for efficient absolute value function + using std::abs; + if((uCoff < 0 && abs(uCoff) > epsilon) || (uCoff > 1 && abs(uCoff - 1.0) > epsilon)) + { + // Ray intersects plane, but not triangle + return {}; + } + + PointT sCrossE1 = sDist.cross(edge1); + + T vCoff = invDet * dirVec.dot(sCrossE1); + + if((vCoff < 0 && abs(vCoff) > epsilon) || (uCoff + vCoff > 1 && abs(uCoff + vCoff - 1.0) > epsilon)) + { + // Ray intersects plane, but not triangle + return {}; + } + + T tCoff = invDet * edge2.dot(sCrossE1); + + if(tCoff > epsilon) + { + // Ray intersection + return {origin + dirVec * tCoff}; + } + + // line intersection not ray intersection + return {}; +} + } // namespace IntersectionUtilities } // namespace nx::core diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 58bd20e886..b53a456762 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -31,6 +31,7 @@ add_executable(simplnx_test GeometryTest.cpp GeometryTestUtilities.hpp H5Test.cpp + IntersectionUtilitiesTest.cpp IOFormat.cpp MontageTest.cpp PluginTest.cpp diff --git a/test/IntersectionUtilitiesTest.cpp b/test/IntersectionUtilitiesTest.cpp new file mode 100644 index 0000000000..e26bb5499d --- /dev/null +++ b/test/IntersectionUtilitiesTest.cpp @@ -0,0 +1,116 @@ +#include "simplnx/Utilities/IntersectionUtilities.hpp" +#include "simplnx/unit_test/simplnx_test_dirs.hpp" + +#include + +#include + +using namespace nx::core; + +TEST_CASE("Simplnx::IntersectionUtilities::Moller-Trumbore Intersection Test (Cached)", "[Simplnx][IntersectionUtilities]") +{ + using PointT = Eigen::Vector3; + IntersectionUtilities::MTPointsCache cache; + cache.pointA = PointT{0.5, 0, 0}; + cache.pointB = PointT{-0.5, 0.5, 0}; + cache.pointC = PointT{-0.5, -0.5, 0}; + + cache.edge1 = cache.pointB - cache.pointA; + cache.edge2 = cache.pointC - cache.pointA; + { + // Case 1 (Intersection) + cache.origin = PointT{0, 0, -1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + REQUIRE(IntersectionUtilities::MTIntersection(PointT{0, 0, 1}, cache).has_value()); + } + + { + // Case 2 (Intersection) + cache.origin = PointT{0, 0, -1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + PointT dirVec = PointT{0.1, 0.1, 1}; + dirVec.normalize(); + + REQUIRE(IntersectionUtilities::MTIntersection(dirVec, cache).has_value()); + } + + { + // Case 3 (No Intersection [Origin Above]) + cache.origin = PointT{0, 0, 1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, 1}, cache).has_value()); + } + + { + // Case 4 (No Intersection [Ray Parallel]) + cache.origin = PointT{0, 0, -1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 1, 0}, cache).has_value()); + } + + { + // Case 5 (No Intersection [Misaligned Origin]) + cache.origin = PointT{10, 0, -1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, 1}, cache).has_value()); + } + + { + // Case 6 (No Intersection [Wrong Direction]) + cache.origin = PointT{0, 0, -1}; + cache.sDist = cache.origin - cache.pointA; + cache.sCrossE1 = cache.sDist.cross(cache.edge1); + + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, -1}, cache).has_value()); + } +} + +TEST_CASE("Simplnx::IntersectionUtilities::Moller-Trumbore Intersection Test", "[Simplnx][IntersectionUtilities]") +{ + using PointT = Eigen::Vector3; + PointT pointA = PointT{0.5, 0, 0}; + PointT pointB = PointT{-0.5, 0.5, 0}; + PointT pointC = PointT{-0.5, -0.5, 0}; + + { + // Case 1 (Intersection) + REQUIRE(IntersectionUtilities::MTIntersection(PointT{0, 0, -1}, PointT{0, 0, 1}, pointA, pointB, pointC).has_value()); + } + + { + // Case 2 (Intersection) + PointT dirVec = PointT{0.1, 0.1, 1}; + dirVec.normalize(); + REQUIRE(IntersectionUtilities::MTIntersection(PointT{0, 0, -1}, dirVec, pointA, pointB, pointC).has_value()); + } + + { + // Case 3 (No Intersection [Origin Above]) + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, 1}, PointT{0, 0, 1}, pointA, pointB, pointC).has_value()); + } + + { + // Case 4 (No Intersection [Ray Parallel]) + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, -1}, PointT{0, 1, 0}, pointA, pointB, pointC).has_value()); + } + + { + // Case 5 (No Intersection [Misaligned Origin]) + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{10, 0, -1}, PointT{0, 0, 1}, pointA, pointB, pointC).has_value()); + } + + { + // Case 6 (No Intersection [Wrong Direction]) + REQUIRE(!IntersectionUtilities::MTIntersection(PointT{0, 0, -1}, PointT{0, 0, -1}, pointA, pointB, pointC).has_value()); + } +}