From 9cf1fb3d80a81c546e0a892faca09703b8791376 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Sun, 1 Jun 2025 20:14:14 -0400 Subject: [PATCH] ENH: Updates codes to use Eigen instead of built in Matrix Operations Signed-off-by: Michael Jackson --- .../Algorithms/GroupMicroTextureRegions.cpp | 91 +++++++++---------- .../Algorithms/GroupMicroTextureRegions.hpp | 8 +- 2 files changed, 45 insertions(+), 54 deletions(-) diff --git a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp index 3ea24c8..7606bde 100644 --- a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp +++ b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp @@ -2,13 +2,14 @@ #include "simplnx/Common/Constants.hpp" #include "simplnx/DataStructure/DataArray.hpp" -#include "simplnx/DataStructure/DataGroup.hpp" #include "simplnx/DataStructure/NeighborList.hpp" #include "simplnx/Utilities/Math/GeometryMath.hpp" #include "simplnx/Utilities/Math/MatrixMath.hpp" #include "EbsdLib/LaueOps/LaueOps.h" +#include "Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/OrientationUtilities.hpp" + #include using namespace nx::core; @@ -19,7 +20,7 @@ void RandomizeFeatureIds(usize totalPoints, usize totalFeatures, Int32Array& cel { // Generate an even distribution of numbers between the min and max range std::mt19937_64 gen(seed); - std::uniform_int_distribution dist(0, totalFeatures - 1); + std::uniform_int_distribution dist(0, static_cast(totalFeatures) - 1); std::vector gid(totalFeatures); std::iota(gid.begin(), gid.end(), 0); @@ -27,7 +28,7 @@ void RandomizeFeatureIds(usize totalPoints, usize totalFeatures, Int32Array& cel //--- Shuffle elements by randomly exchanging each with one other. for(usize i = 1; i < totalFeatures; i++) { - auto r = static_cast(dist(gen)); // Random remaining position. + const auto r = static_cast(dist(gen)); // Random remaining position. if(r >= totalFeatures) { continue; @@ -94,10 +95,10 @@ void GroupMicroTextureRegions::execute() int32 seed = 0; int32 list1size = 0, list2size = 0, listsize = 0; int32 neigh = 0; - bool patchGrouping = false; while(seed >= 0) { + bool patchGrouping = false; parentcount++; seed = getSeed(parentcount); if(seed >= 0) @@ -106,7 +107,7 @@ void GroupMicroTextureRegions::execute() for(std::vector::size_type j = 0; j < grouplist.size(); j++) { int32 firstfeature = grouplist[j]; - list1size = int32(neighborlist[firstfeature].size()); + list1size = static_cast(neighborlist[firstfeature].size()); if(m_InputValues->UseNonContiguousNeighbors) { list2size = nonContigNeighList->getListSize(firstfeature); @@ -156,7 +157,7 @@ void GroupMicroTextureRegions::execute() for(std::vector::size_type j = 0; j < grouplist.size(); j++) { int32_t firstfeature = grouplist[j]; - listsize = int32_t(neighborlist[firstfeature].size()); + listsize = static_cast(neighborlist[firstfeature].size()); for(int32_t l = 0; l < listsize; l++) { neigh = neighborlist[firstfeature][l]; @@ -179,12 +180,10 @@ void GroupMicroTextureRegions::execute() // ----------------------------------------------------------------------------- Result<> GroupMicroTextureRegions::operator()() { + // m_AvgCAxes = {0.0f, 0.0f, 0.0f}; m_Generator = std::mt19937_64(m_InputValues->SeedValue); m_Distribution = std::uniform_real_distribution(0.0f, 1.0f); - m_AvgCAxes[0] = 0.0f; - m_AvgCAxes[1] = 0.0f; - m_AvgCAxes[2] = 0.0f; auto& featureParentIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureParentIdsArrayName); featureParentIds.fill(-1); @@ -205,7 +204,7 @@ Result<> GroupMicroTextureRegions::operator()() cellParentIds[k] = featureParentIds[featureIds[k]]; } - // By default we randomize grains !!! COMMENT OUT FOR DEMONSTRATION !!! + // By default, we randomize grains !!! COMMENT OUT FOR DEMONSTRATION !!! m_MessageHandler(IFilter::Message::Type::Info, "Randomizing Parent Ids"); RandomizeFeatureIds(totalPoints, m_NumTuples, cellParentIds, featureParentIds, featureIds, m_InputValues->SeedValue); @@ -219,9 +218,6 @@ int GroupMicroTextureRegions::getSeed(int32 newFid) auto& featurePhases = m_DataStructure.getDataRefAs(m_InputValues->FeaturePhasesArrayPath); usize numFeatures = featurePhases.getNumberOfTuples(); - - float32 g1[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - float32 g1t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; int32 voxelSeed = -1; // Precalculate some constants @@ -250,6 +246,7 @@ int GroupMicroTextureRegions::getSeed(int32 newFid) // std::ofstream fout ("/tmp/GroupMicroTextureInitialVoxelSeeds.txt", std::ios_base::out | std::ios_base::app); // fout << fmt::format("Feature Parent Id: {} | X: {}, Y: {}\n", voxelSeed, centroids.getComponent(voxelSeed, 0), centroids.getComponent(voxelSeed, 1)); // } + const Eigen::Vector3f cAxis = {0.0f, 0.0f, 1.0f}; if(voxelSeed >= 0) { @@ -261,20 +258,22 @@ int GroupMicroTextureRegions::getSeed(int32 newFid) auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); auto& avgQuats = m_DataStructure.getDataRefAs(m_InputValues->AvgQuatsArrayPath); - usize index = voxelSeed * 4; - OrientationTransformation::qu2om({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g1); - - std::array c1 = {0.0f, 0.0f, 0.0f}; - std::array cAxis = {0.0f, 0.0f, 1.0f}; - // transpose the g matrix so when c-axis is multiplied by it, + const usize index = voxelSeed * 4; + const auto om = OrientationTransformation::qu2om({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}); + // transpose the g matrix so when it multiplies c-axis, // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g1, g1t); - MatrixMath::Multiply3x3with3x1(g1t, cAxis.data(), c1.data()); + const OrientationUtilities::Matrix3fR g1t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om); + + Eigen::Vector3f c1 = g1t * cAxis; // normalize so that the dot product can be taken below without // dividing by the magnitudes (they would be 1) - MatrixMath::Normalize3x1(c1.data()); - MatrixMath::Copy3x1(c1.data(), m_AvgCAxes.data()); - MatrixMath::Multiply3x1withConstant(m_AvgCAxes.data(), volumes.getValue(voxelSeed)); + c1.normalize(); + + c1 = c1 * volumes[voxelSeed]; + + m_RunningAvgOrientation[0] = c1[0]; + m_RunningAvgOrientation[1] = c1[1]; + m_RunningAvgOrientation[2] = c1[2]; } } @@ -284,57 +283,51 @@ int GroupMicroTextureRegions::getSeed(int32 newFid) // ----------------------------------------------------------------------------- bool GroupMicroTextureRegions::determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid) { - uint32 phase1 = 0; - float32 g1[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - float32 g2[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - float32 g1t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - float32 g2t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}; - std::array c1 = {0.0f, 0.0f, 0.0f}; - std::array caxis = {0.0f, 0.0f, 1.0f}; + Eigen::Vector3f cAxis = {0.0f, 0.0f, 1.0f}; auto& featurePhases = m_DataStructure.getDataRefAs(m_InputValues->FeaturePhasesArrayPath); auto& featureParentIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureParentIdsArrayName); auto& crystalStructures = m_DataStructure.getDataRefAs(m_InputValues->CrystalStructuresArrayPath); if(featureParentIds[neighborFeature] == -1 && featurePhases[referenceFeature] > 0 && featurePhases[neighborFeature] > 0) { + uint32 phase1 = 0; + Eigen::Vector3f c1; auto& avgQuats = m_DataStructure.getDataRefAs(m_InputValues->AvgQuatsArrayPath); if(!m_InputValues->UseRunningAverage) { usize index = referenceFeature * 4; phase1 = crystalStructures[featurePhases[referenceFeature]]; - OrientationTransformation::qu2om>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g1); - - // transpose the g matrix so when c-axis is multiplied by it, + const auto om = OrientationTransformation::qu2om({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}); + // transpose the g matrix so when it multiplies c-axis, // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g1, g1t); - MatrixMath::Multiply3x3with3x1(g1t, caxis.data(), c1.data()); + const OrientationUtilities::Matrix3fR g1t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om); + c1 = g1t * cAxis; // normalize so that the dot product can be taken below without // dividing by the magnitudes (they would be 1) - MatrixMath::Normalize3x1(c1.data()); + c1.normalize(); } uint32 phase2 = crystalStructures[featurePhases[neighborFeature]]; if(phase1 == phase2 && (phase1 == EbsdLib::CrystalStructure::Hexagonal_High)) { usize index = neighborFeature * 4; - OrientationTransformation::qu2om({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g2); - - std::array c2 = {0.0f, 0.0f, 0.0f}; - // transpose the g matrix so when c-axis is multiplied by it, + const auto om = OrientationTransformation::qu2om({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}); + // transpose the g matrix so when it multiplies c-axis, // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g2, g2t); - MatrixMath::Multiply3x3with3x1(g2t, caxis.data(), c2.data()); + const OrientationUtilities::Matrix3fR g2t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om); + + Eigen::Vector3f c2 = g2t * cAxis; // normalize so that the dot product can be taken below without // dividing by the magnitudes (they would be 1) - MatrixMath::Normalize3x1(c2.data()); + c2.normalize(); float32 w; if(m_InputValues->UseRunningAverage) { - w = GeometryMath::CosThetaBetweenVectors(Point3Df{m_AvgCAxes}, Point3Df{c2}); + w = GeometryMath::CosThetaBetweenVectors(m_RunningAvgOrientation, c2); } else { - w = GeometryMath::CosThetaBetweenVectors(Point3Df{c1}, Point3Df{c2}); + w = GeometryMath::CosThetaBetweenVectors(c1, c2); } if(w < -1.0f) @@ -348,15 +341,15 @@ bool GroupMicroTextureRegions::determineGrouping(int32 referenceFeature, int32 n w = std::acos(w); // Convert user defined tolerance to radians. - float32 cAxisToleranceRad = m_InputValues->CAxisTolerance * nx::core::Constants::k_PiD / 180.0f; + float32 cAxisToleranceRad = m_InputValues->CAxisTolerance * nx::core::Constants::k_PiF / 180.0f; if(w <= cAxisToleranceRad || (nx::core::Constants::k_PiD - w) <= cAxisToleranceRad) { featureParentIds[neighborFeature] = newFid; if(m_InputValues->UseRunningAverage) { auto& volumes = m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath); - MatrixMath::Multiply3x1withConstant(c2.data(), volumes.getValue(neighborFeature)); - MatrixMath::Add3x1s(m_AvgCAxes.data(), c2.data(), m_AvgCAxes.data()); + MatrixMath::Multiply3x1withConstant(c2.data(), volumes[neighborFeature]); + MatrixMath::Add3x1s(m_RunningAvgOrientation.data(), c2.data(), m_RunningAvgOrientation.data()); } return true; } diff --git a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp index 247116c..7367f1d 100644 --- a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp +++ b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp @@ -5,10 +5,8 @@ #include "simplnx/DataStructure/DataPath.hpp" #include "simplnx/DataStructure/DataStructure.hpp" #include "simplnx/Filter/IFilter.hpp" -#include "simplnx/Parameters/ArraySelectionParameter.hpp" -#include "simplnx/Parameters/BoolParameter.hpp" -#include "simplnx/Parameters/NumberParameter.hpp" -#include "simplnx/Parameters/StringParameter.hpp" + +#include #include @@ -67,7 +65,7 @@ class SIMPLNXREVIEW_EXPORT GroupMicroTextureRegions const IFilter::MessageHandler& m_MessageHandler; usize m_NumTuples = 0; - std::array m_AvgCAxes = {0.0f, 0.0f, 0.0f}; + Eigen::Vector3 m_RunningAvgOrientation = {0.0f, 0.0f, 0.0f}; std::mt19937_64 m_Generator = {}; std::uniform_real_distribution m_Distribution = {}; };