ENH: Fixes algorithm to compute shape factors from Triangle geometry (B…

…lueQuartzSoftware#1157)

* Add documentation to the functions.

All values are validated against the voxel generated results and the theoretical values.

Theoretical values taken from paper from De Graef, Simmons, et. al.

Signed-off-by: Michael Jackson <[email protected]>

* Use eigenvalues and eigenvector to find principle axes and primary axes

* correct production of eigenvalues/vectors

* Correct axis length production, spacial requirement reduction, code consolidation, further comment documentation

* parameter update: remove unused volumes array

* Patch test and calculate Euler Characteristic fast watertight validation

* use unordered set in fast edge counter, and move to utility file for shared use

* Implement MT algorithm

* Intersection Utilities updated with MT algorithm Cached/Uncached, added test cases for them

* Filter cleanup and optimization

* documentation update

* vastly improved error message

* Test and clang formatting

---------

Signed-off-by: Michael Jackson <[email protected]>
Co-authored-by: nyoungbq <[email protected]>
Co-authored-by: Jared Duffey <[email protected]>

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

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

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.

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<ImageGeom>(m_InputValues->ImageGeometryPath);
 
   const auto& featureIds = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeatureIdsArrayPath);
   const auto& centroids = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->CentroidsArrayPath);
+  // Calculated Arrays
   auto& volumes = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->VolumesArrayPath);
   auto& omega3s = m_DataStructure.getDataRefAs<Float32Array>(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<float>(m_ScaleFactor)));
         ydist1 = (y1 - (centroids[gnum * 3 + 1] * static_cast<float>(m_ScaleFactor)));
         zdist1 = (z1 - (centroids[gnum * 3 + 2] * static_cast<float>(m_ScaleFactor)));
+
         xdist2 = (x1 - (centroids[gnum * 3 + 0] * static_cast<float>(m_ScaleFactor)));
         ydist2 = (y1 - (centroids[gnum * 3 + 1] * static_cast<float>(m_ScaleFactor)));
         zdist2 = (z2 - (centroids[gnum * 3 + 2] * static_cast<float>(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<float>(-m_FeatureMoments[featureId * 6 + 3]);
     u011 = static_cast<float>(-m_FeatureMoments[featureId * 6 + 4]);
     u101 = static_cast<float>(-m_FeatureMoments[featureId * 6 + 5]);
+
     o3 = static_cast<double>((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<Int32Array>(m_InputValues->FeatureIdsArrayPath);
@@ -430,9 +433,7 @@ void ComputeShapes::find_moments2D()
 }
 
 // -----------------------------------------------------------------------------
-//
-// -----------------------------------------------------------------------------
-void ComputeShapes::find_axes()
+void ComputeShapes::findAxes()
 {
   auto& axisLengths = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AxisLengthsArrayPath);
   auto& aspectRatios = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AspectRatiosArrayPath);
@@ -476,9 +477,7 @@ void ComputeShapes::find_axes()
 }
 
 // -----------------------------------------------------------------------------
-//
-// -----------------------------------------------------------------------------
-void ComputeShapes::find_axes2D()
+void ComputeShapes::findAxes2D()
 {
   auto& volumes = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->VolumesArrayPath);
   auto& axisLengths = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AxisLengthsArrayPath);
@@ -551,9 +550,7 @@ void ComputeShapes::find_axes2D()
 }
 
 // -----------------------------------------------------------------------------
-//
-// -----------------------------------------------------------------------------
-void ComputeShapes::find_axiseulers()
+void ComputeShapes::findAxisEulers()
 {
   const auto& centroids = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->CentroidsArrayPath);
   auto& axisEulerAngles = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AxisEulerAnglesArrayPath);
@@ -588,9 +585,7 @@ void ComputeShapes::find_axiseulers()
 }
 
 // -----------------------------------------------------------------------------
-//
-// -----------------------------------------------------------------------------
-void ComputeShapes::find_axiseulers2D()
+void ComputeShapes::findAxisEulers2D()
 {
   const auto& centroids = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->CentroidsArrayPath);
   auto& axisEulerAngles = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AxisEulerAnglesArrayPath);

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