diff --git a/docs/MergeColoniesFilter.md b/docs/MergeColoniesFilter.md index a14c913..891d7ea 100644 --- a/docs/MergeColoniesFilter.md +++ b/docs/MergeColoniesFilter.md @@ -9,23 +9,33 @@ Reconstruction (Grouping) ## Description ## -This **Filter** groups neighboring **Features** that have a *special* misorientation that is associated with *alpha* variants that transformed from the same *beta* grain in titanium. The algorithm for grouping the **Features** is analogous to the algorithm for segmenting the **Features**, except the average orientation of the **Features** are used instead of the orientations of the individual **Elements** and the criterion for grouping is specific to the *alpha-beta transformation*. The user can specify a tolerance on both the *axis* and the *angle* that defines the misorientation relationship (i.e., a tolerance of 1 degree for both tolerances would allow the neighboring **Features** to be grouped if their misorientation was between 59-61 degrees about an axis within 1 degree of a2, as given by the 3rd *special* misorientation below). +This **Filter** groups neighboring **Features** that have a *special* misorientation that is associated with +*alpha* variants that transformed from the same *beta* grain in titanium. The algorithm for grouping the +**Features** is analogous to the algorithm for segmenting the **Features**, except the average orientation +of the **Features** are used instead of the orientations of the individual **Elements** and the criterion +for grouping is specific to the *alpha-beta transformation*. The user can specify a tolerance on both the +*axis* and the *angle* that defines the misorientation relationship (i.e., a tolerance of 1 degree for both +tolerances would allow the neighboring **Features** to be grouped if their misorientation was between 59-61 +degrees about an axis within 1 degree of a2, as given by the 3rd *special* misorientation below). The list of *special* misorientations can be found in the paper by Germain et al.1 and are listed here: -| Angle | Axis | -|------|------| -| 0 | Identity | -| 10.529 | c = <0001> | -| 60 | a2 = <-12-10> | +| Angle | Axis | +|--------|---------------------------------------------------| +| 0 | Identity | +| 10.529 | c = <0001> | +| 60 | a2 = <-12-10> | | 60.832 | d1 at 80.97 degrees from c in the plane of (d3,c) | | 63.262 | d2 at 72.73 degrees from c in the plane of (a2,c) | -| 90 | d3 at 5.26 degrees from a2 in the basal plane | +| 90 | d3 at 5.26 degrees from a2 in the basal plane | % Auto generated parameter table will be inserted here ## References +[1] L. Germain, N. Gey and M. Humbert, Reliability of reconstructed Beta orientation maps in titanium alloys, Ultramicrscopy, 2007, 1129-1135. + + ## Example Pipelines ## License & Copyright diff --git a/src/SimplnxReview/Filters/Algorithms/ComputeGroupingDensity.cpp b/src/SimplnxReview/Filters/Algorithms/ComputeGroupingDensity.cpp index 5d80379..31e20a2 100644 --- a/src/SimplnxReview/Filters/Algorithms/ComputeGroupingDensity.cpp +++ b/src/SimplnxReview/Filters/Algorithms/ComputeGroupingDensity.cpp @@ -2,6 +2,7 @@ #include "simplnx/DataStructure/DataArray.hpp" #include "simplnx/DataStructure/NeighborList.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" using namespace nx::core; @@ -66,30 +67,13 @@ class FindDensityGrouping // Default value-initialized to zeroes: https://en.cppreference.com/w/cpp/named_req/DefaultInsertable checkedFeatureVolumes.resize(numFeatures); } - int32_t progInt = 0; - usize prevParentId = 1; - usize currentParentId = 1; - auto start = std::chrono::steady_clock::now(); + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); for(usize parentIdx = 1; parentIdx < numParents; parentIdx++) { - progInt = static_cast(parentIdx) / static_cast(numParents) * 100.0f; - auto now = std::chrono::steady_clock::now(); - // Only send updates every 1 second - if(std::chrono::duration_cast(now - start).count() > 1000) - { - currentParentId = parentIdx; - auto totalParentIds = currentParentId - prevParentId; - auto rate = static_cast(totalParentIds) / static_cast(std::chrono::duration_cast(now - start).count()); - - auto remainingParents = numParents - parentIdx; - auto minutesRemain = (remainingParents / rate) / 60; // Convert to minutes + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("[{}%]", CalculatePercentComplete(parentIdx, numParents)); }); - std::string message = fmt::format("{}/{} [{}%] at {} parents/sec. Time Remain: {:.2f} Minutes", parentIdx, numParents, progInt, rate, minutesRemain); - m_MessageHandler(nx::core::IFilter::ProgressMessage{nx::core::IFilter::Message::Type::Info, message, progInt}); - start = std::chrono::steady_clock::now(); - prevParentId = currentParentId; - } if(m_ShouldCancel) { return {}; diff --git a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp index 36a0187..8e5e688 100644 --- a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp +++ b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp @@ -2,10 +2,10 @@ #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 "simplnx/Utilities/MessageHelper.hpp" #include "EbsdLib/LaueOps/LaueOps.h" @@ -13,40 +13,6 @@ using namespace nx::core; -namespace -{ -void RandomizeFeatureIds(usize totalPoints, usize totalFeatures, Int32Array& cellParentIds, Int32Array& featureParentIds, const Int32Array& featureIds, uint64 seed) -{ - // 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::vector gid(totalFeatures); - std::iota(gid.begin(), gid.end(), 0); - - //--- 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. - if(r >= totalFeatures) - { - continue; - } - - int32 temp = gid[i]; - gid[i] = gid[r]; - gid[r] = temp; - } - - // Now adjust all the Grain id values for each Voxel - for(usize i = 0; i < totalPoints; ++i) - { - cellParentIds[i] = gid[cellParentIds[i]]; - featureParentIds[featureIds[i]] = cellParentIds[i]; - } -} -} // namespace - // ----------------------------------------------------------------------------- GroupMicroTextureRegions::GroupMicroTextureRegions(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, GroupMicroTextureRegionsInputValues* inputValues) @@ -54,6 +20,11 @@ GroupMicroTextureRegions::GroupMicroTextureRegions(DataStructure& dataStructure, , m_InputValues(inputValues) , m_ShouldCancel(shouldCancel) , m_MessageHandler(mesgHandler) +, m_FeaturePhases(m_DataStructure.getDataRefAs(m_InputValues->FeaturePhasesArrayPath)) +, m_FeatureParentIds(m_DataStructure.getDataRefAs(m_InputValues->FeatureParentIdsArrayName)) +, m_CrystalStructures(m_DataStructure.getDataRefAs(m_InputValues->CrystalStructuresArrayPath)) +, m_AvgQuats(m_DataStructure.getDataRefAs(m_InputValues->AvgQuatsArrayPath)) +, m_Volumes(m_DataStructure.getDataRefAs(m_InputValues->VolumesArrayPath)) { } @@ -79,37 +50,44 @@ bool GroupMicroTextureRegions::growGrouping(int32_t referenceFeature, int32_t ne } // ----------------------------------------------------------------------------- -void GroupMicroTextureRegions::execute() +Result<> GroupMicroTextureRegions::execute() { - NeighborList& neighborlist = m_DataStructure.getDataRefAs>(m_InputValues->ContiguousNeighborListArrayPath); - NeighborList* nonContigNeighList = nullptr; + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + + NeighborList& featureNeighborListRef = m_DataStructure.getDataRefAs>(m_InputValues->ContiguousNeighborListArrayPath); + NeighborList* nonContigNeighListPtr = nullptr; if(m_InputValues->UseNonContiguousNeighbors) { - nonContigNeighList = m_DataStructure.getDataAs>(m_InputValues->NonContiguousNeighborListArrayPath); + nonContigNeighListPtr = m_DataStructure.getDataAs>(m_InputValues->NonContiguousNeighborListArrayPath); + } + if(nullptr == nonContigNeighListPtr) + { + return MakeErrorResult(-99345, "There was an error getting the Non-contiguous neighborlist from the DataStructure"); } - std::vector grouplist; + std::vector groupList; + + int32 parentCount = 0; + int32 featureSeed = 0; + int32 list1size = 0, list2size = 0, listSize = 0; - int32 parentcount = 0; - int32 seed = 0; - int32 list1size = 0, list2size = 0, listsize = 0; - int32 neigh = 0; bool patchGrouping = false; - while(seed >= 0) + while(featureSeed >= 0) { - parentcount++; - seed = getSeed(parentcount); - if(seed >= 0) + parentCount++; + featureSeed = getSeed(parentCount); + if(featureSeed >= 0) { - grouplist.push_back(seed); - for(std::vector::size_type j = 0; j < grouplist.size(); j++) + groupList.push_back(featureSeed); + for(std::vector::size_type j = 0; j < groupList.size(); j++) { - int32 firstfeature = grouplist[j]; - list1size = int32(neighborlist[firstfeature].size()); + int32 firstFeature = groupList[j]; + list1size = static_cast(featureNeighborListRef[firstFeature].size()); if(m_InputValues->UseNonContiguousNeighbors) { - list2size = nonContigNeighList->getListSize(firstfeature); + list2size = nonContigNeighListPtr->getListSize(firstFeature); } for(int32 k = 0; k < 2; k++) { @@ -119,30 +97,31 @@ void GroupMicroTextureRegions::execute() } if(k == 0) { - listsize = list1size; + listSize = list1size; } else if(k == 1) { - listsize = list2size; + listSize = list2size; } - for(int32 l = 0; l < listsize; l++) + for(int32 l = 0; l < listSize; l++) { + int32 neigh = -1; if(k == 0) { - neigh = neighborlist[firstfeature][l]; + neigh = featureNeighborListRef[firstFeature][l]; } - else if(k == 1) + else if(k == 1 && m_InputValues->UseNonContiguousNeighbors) { bool ok = false; - neigh = nonContigNeighList->getValue(firstfeature, l, ok); + neigh = nonContigNeighListPtr->getValue(firstFeature, l, ok); } - if(neigh != firstfeature) + if(neigh >= 0 && neigh != firstFeature) { - if(determineGrouping(firstfeature, neigh, parentcount)) + if(determineGrouping(firstFeature, neigh, parentCount)) { if(!patchGrouping) { - grouplist.push_back(neigh); + groupList.push_back(neigh); } } } @@ -151,44 +130,57 @@ void GroupMicroTextureRegions::execute() } if(patchGrouping) { - if(growPatch(parentcount)) + if(growPatch(parentCount)) { - for(std::vector::size_type j = 0; j < grouplist.size(); j++) + for(std::vector::size_type j = 0; j < groupList.size(); j++) { - int32_t firstfeature = grouplist[j]; - listsize = int32_t(neighborlist[firstfeature].size()); - for(int32_t l = 0; l < listsize; l++) + int32_t firstFeature = groupList[j]; + listSize = static_cast(featureNeighborListRef[firstFeature].size()); + for(int32_t l = 0; l < listSize; l++) { - neigh = neighborlist[firstfeature][l]; - if(neigh != firstfeature) + int32 neigh = featureNeighborListRef[firstFeature][l]; + if(neigh != firstFeature) { - if(growGrouping(firstfeature, neigh, parentcount)) + if(growGrouping(firstFeature, neigh, parentCount)) { - grouplist.push_back(neigh); + groupList.push_back(neigh); } } } } } } + + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Parent Count: {}", parentCount); }); } - grouplist.clear(); + groupList.clear(); } + return {}; } // ----------------------------------------------------------------------------- Result<> GroupMicroTextureRegions::operator()() { - m_Generator = std::mt19937_64(m_InputValues->SeedValue); + MessageHelper messageHelper(m_MessageHandler); + + m_Generator = std::mt19937_64(std::mt19937::default_seed); m_Distribution = std::uniform_real_distribution(0.0f, 1.0f); + // Initialize Data 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); + m_FeatureParentIds.fill(-1); - execute(); + // Execute the main grouping algorithm + messageHelper.sendMessage(fmt::format("Starting Grouping.....")); + + // Execute the grouping algorithm + Result<> result = execute(); + if(result.invalid()) + { + return result; + } // handle active array resize if(m_NumTuples < 2) @@ -202,12 +194,12 @@ Result<> GroupMicroTextureRegions::operator()() usize totalPoints = featureIds.getNumberOfTuples(); for(usize k = 0; k < totalPoints; k++) { - cellParentIds[k] = featureParentIds[featureIds[k]]; + cellParentIds[k] = m_FeatureParentIds[featureIds[k]]; } - // 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); + // By default, we randomize grains !!! COMMENT OUT FOR DEMONSTRATION !!! + // m_MessageHandler(IFilter::Message::Type::Info, "Randomizing Parent Ids"); + // RandomizeFeatureIds(totalPoints, m_NumTuples, cellParentIds, m_FeatureParentIds, featureIds, m_InputValues->SeedValue); return {}; } @@ -215,29 +207,32 @@ Result<> GroupMicroTextureRegions::operator()() // ----------------------------------------------------------------------------- int GroupMicroTextureRegions::getSeed(int32 newFid) { - auto& featureParentIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureParentIdsArrayName); - auto& featurePhases = m_DataStructure.getDataRefAs(m_InputValues->FeaturePhasesArrayPath); - - usize numFeatures = featurePhases.getNumberOfTuples(); + usize numFeatures = m_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; + int32 featureIdSeed = -1; // Precalculate some constants - int32 totalFMinus1 = numFeatures - 1; + const int32 totalFMinus1 = static_cast(numFeatures) - 1; usize counter = 0; + // This section finds a feature id that has not been grouped yet. It starts by + // randomly selecting a feature id between 0 and numFeatures-1. We then start + // looping. If the initial random value is valid then we exit the loop after + // a single iteration. If that feature has already been grouped, then we add one + // to the `randFeature` value and try again. If we get to the end of the range of + // featureIds then the algorithm will loop back to featureId = 0 and start incrementing + // from there. This is reasonably efficient as we only generate random numbers + // as needed. auto randFeature = static_cast(m_Distribution(m_Generator) * static_cast(totalFMinus1)); - while(voxelSeed == -1 && counter < numFeatures) + while(featureIdSeed == -1 && counter < numFeatures) { if(randFeature > totalFMinus1) { randFeature = randFeature - numFeatures; } - if(featureParentIds[randFeature] == -1) + if(m_FeatureParentIds.getValue(randFeature) == -1) { - voxelSeed = randFeature; + featureIdSeed = randFeature; } randFeature++; counter++; @@ -251,112 +246,89 @@ int GroupMicroTextureRegions::getSeed(int32 newFid) // fout << fmt::format("Feature Parent Id: {} | X: {}, Y: {}\n", voxelSeed, centroids.getComponent(voxelSeed, 0), centroids.getComponent(voxelSeed, 1)); // } - if(voxelSeed >= 0) + if(featureIdSeed >= 0) { - featureParentIds[voxelSeed] = newFid; + m_FeatureParentIds[featureIdSeed] = newFid; m_NumTuples = newFid + 1; if(m_InputValues->UseRunningAverage) { - 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, - // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g1, g1t); - MatrixMath::Multiply3x3with3x1(g1t, cAxis.data(), c1.data()); + usize index = featureIdSeed * 4; + // Get the orientation matrix (which is passive) and then transpose it to make it active transform + ebsdlib::Matrix3X3F g1t = ebsdlib::Quaternion(m_AvgQuats.getValue(index + 0), m_AvgQuats.getValue(index + 1), m_AvgQuats.getValue(index + 2), m_AvgQuats.getValue(index + 3)) + .toOrientationMatrix() + .toGMatrix() + .transpose(); + ebsdlib::Matrix3X1F cAxis(0.0f, 0.0f, 1.0f); // 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)); + const ebsdlib::Matrix3X1F c1 = (g1t * cAxis).normalize(); + + m_AvgCAxes = c1 * m_Volumes.getValue(featureIdSeed); } } - return voxelSeed; + return featureIdSeed; } // ----------------------------------------------------------------------------- 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}; - - 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) + const int32 neighborParentId = m_FeatureParentIds.getValue(neighborFeature); + const int32 referenceFeaturePhase = m_FeaturePhases.getValue(referenceFeature); + const int32 neighborFeaturePhase = m_FeaturePhases.getValue(neighborFeature); + + if(neighborParentId == -1 && referenceFeaturePhase > 0 && neighborFeaturePhase > 0) { - auto& avgQuats = m_DataStructure.getDataRefAs(m_InputValues->AvgQuatsArrayPath); + ebsdlib::Matrix3X1F c1 = {0.0f, 0.0f, 0.0f}; + ebsdlib::Matrix3X1F cAxis(0.0f, 0.0f, 1.0f); + 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); - + const usize index = referenceFeature * 4; + // Get the orientation matrix (which is passive) and then transpose it to make it active transform // transpose the g matrix so when c-axis is multiplied by it, // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g1, g1t); - MatrixMath::Multiply3x3with3x1(g1t, caxis.data(), c1.data()); - // normalize so that the dot product can be taken below without - // dividing by the magnitudes (they would be 1) - MatrixMath::Normalize3x1(c1.data()); + ebsdlib::Matrix3X3F g1t = ebsdlib::Quaternion(m_AvgQuats.getValue(index + 0), m_AvgQuats.getValue(index + 1), m_AvgQuats.getValue(index + 2), m_AvgQuats.getValue(index + 3)) + .toOrientationMatrix() + .toGMatrix() + .transpose(); + c1 = (g1t * cAxis).normalize(); } - uint32 phase2 = crystalStructures[featurePhases[neighborFeature]]; - if(phase1 == phase2 && (phase1 == EbsdLib::CrystalStructure::Hexagonal_High)) + uint32 phase1 = m_CrystalStructures.getValue(referenceFeaturePhase); + uint32 phase2 = m_CrystalStructures.getValue(neighborFeaturePhase); + 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}; + const usize index = neighborFeature * 4; + // Get the orientation matrix (which is passive) and then transpose it to make it active transform // transpose the g matrix so when c-axis is multiplied by it, // it will give the sample direction that the c-axis is along - MatrixMath::Transpose3x3(g2, g2t); - MatrixMath::Multiply3x3with3x1(g2t, caxis.data(), c2.data()); - // normalize so that the dot product can be taken below without - // dividing by the magnitudes (they would be 1) - MatrixMath::Normalize3x1(c2.data()); + ebsdlib::Matrix3X3F g2t = ebsdlib::Quaternion(m_AvgQuats.getValue(index + 0), m_AvgQuats.getValue(index + 1), m_AvgQuats.getValue(index + 2), m_AvgQuats.getValue(index + 3)) + .toOrientationMatrix() + .toGMatrix() + .transpose(); + ebsdlib::Matrix3X1F c2 = (g2t * cAxis).normalize(); float32 w; if(m_InputValues->UseRunningAverage) { - w = GeometryMath::CosThetaBetweenVectors(Point3Df{m_AvgCAxes}, Point3Df{c2}); + w = m_AvgCAxes.cosTheta(c2); } else { - w = GeometryMath::CosThetaBetweenVectors(Point3Df{c1}, Point3Df{c2}); - } - - if(w < -1.0f) - { - w = -1.0f; - } - else if(w > 1.0f) - { - w = 1.0f; + w = c1.cosTheta(c2); } - w = std::acos(w); + w = std::acos(std::clamp(w, -1.0f, 1.0f)); // 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; + m_FeatureParentIds.setValue(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()); + c2 = c2 * m_Volumes.getValue(neighborFeature); + m_AvgCAxes = m_AvgCAxes + c2; } return true; } diff --git a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp index 247116c..923924f 100644 --- a/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp +++ b/src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp @@ -2,13 +2,12 @@ #include "SimplnxReview/SimplnxReview_export.hpp" +#include "simplnx/DataStructure/DataArray.hpp" #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 "EbsdLib/Math/Matrix3X1.hpp" #include @@ -56,7 +55,7 @@ class SIMPLNXREVIEW_EXPORT GroupMicroTextureRegions protected: int getSeed(int32 newFid); bool determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid); - void execute(); + Result<> execute(); bool growPatch(int32 currentPatch); bool growGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid); @@ -67,8 +66,16 @@ class SIMPLNXREVIEW_EXPORT GroupMicroTextureRegions const IFilter::MessageHandler& m_MessageHandler; usize m_NumTuples = 0; - std::array m_AvgCAxes = {0.0f, 0.0f, 0.0f}; + ebsdlib::Matrix3X1F m_AvgCAxes = {0.0f, 0.0f, 0.0f}; std::mt19937_64 m_Generator = {}; std::uniform_real_distribution m_Distribution = {}; + + // These are so that we don't have to keep getting the references while we are running + + Int32Array& m_FeaturePhases; + Int32Array& m_FeatureParentIds; + UInt32Array& m_CrystalStructures; + Float32Array& m_AvgQuats; + Float32Array& m_Volumes; }; } // namespace nx::core diff --git a/src/SimplnxReview/Filters/Algorithms/MergeColonies.cpp b/src/SimplnxReview/Filters/Algorithms/MergeColonies.cpp index a07660c..11657f1 100644 --- a/src/SimplnxReview/Filters/Algorithms/MergeColonies.cpp +++ b/src/SimplnxReview/Filters/Algorithms/MergeColonies.cpp @@ -3,18 +3,18 @@ #include "simplnx/Common/Array.hpp" #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 "simplnx/Utilities/MessageHelper.hpp" #include "EbsdLib/Core/EbsdLibConstants.h" #include "EbsdLib/Core/Orientation.hpp" -#include "EbsdLib/Core/OrientationTransformation.hpp" -#include "EbsdLib/Core/Quaternion.hpp" +#include "EbsdLib/Orientation/Quaternion.hpp" + +#include using namespace nx::core; -using LaueOpsShPtrType = std::shared_ptr; +using LaueOpsShPtrType = std::shared_ptr; using LaueOpsContainer = std::vector; namespace @@ -24,70 +24,51 @@ const float64 unit111 = 1.0 / std::sqrt(3.0); const float64 unit112_1 = 1.0 / std::sqrt(6.0); const float64 unit112_2 = 2.0 / std::sqrt(6.0); -float64 crystalDirections[12][3][3] = {{{unit111, unit112_1, unit110}, {-unit111, -unit112_1, unit110}, {unit111, -unit112_2, 0}}, - - {{-unit111, unit112_1, unit110}, {unit111, -unit112_1, unit110}, {unit111, unit112_2, 0}}, - - {{unit111, -unit112_1, unit110}, {unit111, -unit112_1, -unit110}, {unit111, unit112_2, 0}}, - - {{unit111, unit112_1, unit110}, {unit111, unit112_1, -unit110}, {-unit111, unit112_2, 0}}, - - {{unit111, unit112_1, unit110}, {unit111, -unit112_2, 0}, {unit111, unit112_1, -unit110}}, - - {{unit111, -unit112_1, unit110}, {-unit111, -unit112_2, 0}, {unit111, -unit112_1, -unit110}}, - - {{unit111, -unit112_1, unit110}, {unit111, unit112_2, 0}, {-unit111, unit112_1, unit110}}, - - {{-unit111, -unit112_1, unit110}, {unit111, -unit112_2, 0}, {unit111, unit112_1, unit110}}, - - {{unit111, -unit112_2, 0}, {unit111, unit112_1, unit110}, {-unit111, -unit112_1, unit110}}, - - {{unit111, unit112_2, 0}, {-unit111, unit112_1, unit110}, {unit111, -unit112_1, unit110}}, - - {{unit111, unit112_2, 0}, {unit111, -unit112_1, unit110}, {unit111, -unit112_1, -unit110}}, - - {{-unit111, unit112_2, 0}, {unit111, unit112_1, unit110}, {unit111, unit112_1, -unit110}}}; +std::vector crystalDirections = { + {unit111, unit112_1, unit110, -unit111, -unit112_1, unit110, unit111, -unit112_2, 0}, {-unit111, unit112_1, unit110, unit111, -unit112_1, unit110, unit111, unit112_2, 0}, + {unit111, -unit112_1, unit110, unit111, -unit112_1, -unit110, unit111, unit112_2, 0}, {unit111, unit112_1, unit110, unit111, unit112_1, -unit110, -unit111, unit112_2, 0}, + {unit111, unit112_1, unit110, unit111, -unit112_2, 0, unit111, unit112_1, -unit110}, {unit111, -unit112_1, unit110, -unit111, -unit112_2, 0, unit111, -unit112_1, -unit110}, + {unit111, -unit112_1, unit110, unit111, unit112_2, 0, -unit111, unit112_1, unit110}, {-unit111, -unit112_1, unit110, unit111, -unit112_2, 0, unit111, unit112_1, unit110}, + {unit111, -unit112_2, 0, unit111, unit112_1, unit110, -unit111, -unit112_1, unit110}, {unit111, unit112_2, 0, -unit111, unit112_1, unit110, unit111, -unit112_1, unit110}, + {unit111, unit112_2, 0, unit111, -unit112_1, unit110, unit111, -unit112_1, -unit110}, {-unit111, unit112_2, 0, unit111, unit112_1, unit110, unit111, unit112_1, -unit110}}; // ----------------------------------------------------------------------------- // // ----------------------------------------------------------------------------- -bool check_for_burgers(const QuatD& betaQuat, const QuatD& alphaQuat, float64 angleTolerance) +bool check_for_burgers(const ebsdlib::QuatD& betaQuat, const ebsdlib::QuatD& alphaQuat, float64 angleTolerance) { float64 dP = 0.0; float64 angle = 0.0; - float64 radToDeg = 180.0 / Constants::k_PiD; + constexpr float64 radToDeg = 180.0 / Constants::k_PiD; - float64 gBeta[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - float64 gBetaT[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - OrientationTransformation::qu2om(betaQuat).toGMatrix(gBeta); // transpose gBeta so the sample direction is the output when // gBeta is multiplied by the crystal directions below - MatrixMath::Transpose3x3(gBeta, gBetaT); - - float64 gAlpha[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - float64 gAlphaT[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; - OrientationTransformation::qu2om(alphaQuat).toGMatrix(gAlpha); + const ebsdlib::Matrix3X3D gBetaT = betaQuat.toOrientationMatrix().transpose().toGMatrix(); // transpose gBeta so the sample direction is the output when // gBeta is multiplied by the crystal directions below - MatrixMath::Transpose3x3(gAlpha, gAlphaT); + ebsdlib::Matrix3X3D gAlphaT = alphaQuat.toOrientationMatrix().transpose().toGMatrix(); - float64 mat[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; for(int32 i = 0; i < 12; i++) { - MatrixMath::Multiply3x3with3x3(gBetaT, crystalDirections[i], mat); - Point3Dd a = Point3Dd(mat[0][2], mat[1][2], mat[2][2]); - Point3Dd b = Point3Dd(gAlphaT[0][2], gAlphaT[1][2], gAlphaT[2][2]); - dP = GeometryMath::CosThetaBetweenVectors(a, b); + ebsdlib::Matrix3X3D mat = gBetaT * crystalDirections[i]; + + ebsdlib::Matrix3X1D a(mat[2], mat[5], mat[8]); + ebsdlib::Matrix3X1D b(gAlphaT[2], gAlphaT[5], gAlphaT[8]); + + dP = a.cosTheta(b); + dP = std::clamp(dP, -1.0, 1.0); angle = std::acos(dP); + if((angle * radToDeg) < angleTolerance || (180.0f - (angle * radToDeg)) < angleTolerance) { - a[0] = mat[0][0]; - a[1] = mat[1][0]; - a[2] = mat[2][0]; - b[0] = gAlphaT[0][0]; - b[1] = gAlphaT[1][0]; - b[2] = gAlphaT[2][0]; - dP = GeometryMath::CosThetaBetweenVectors(a, b); + a[0] = mat[0]; + a[1] = mat[3]; + a[2] = mat[6]; + b[0] = gAlphaT[0]; + b[1] = gAlphaT[3]; + b[2] = gAlphaT[6]; + dP = a.cosTheta(b); + dP = std::clamp(dP, -1.0, 1.0); angle = std::acos(dP); if((angle * radToDeg) < angleTolerance) { @@ -97,10 +78,11 @@ bool check_for_burgers(const QuatD& betaQuat, const QuatD& alphaQuat, float64 an { return true; } - b[0] = -0.5 * gAlphaT[0][0] + 0.866025 * gAlphaT[0][1]; - b[1] = -0.5 * gAlphaT[1][0] + 0.866025 * gAlphaT[1][1]; - b[2] = -0.5 * gAlphaT[2][0] + 0.866025 * gAlphaT[2][1]; - dP = GeometryMath::CosThetaBetweenVectors(a, b); + b[0] = -0.5 * gAlphaT[0] + 0.866025 * gAlphaT[1]; + b[1] = -0.5 * gAlphaT[3] + 0.866025 * gAlphaT[4]; + b[2] = -0.5 * gAlphaT[6] + 0.866025 * gAlphaT[7]; + dP = a.cosTheta(b); + dP = std::clamp(dP, -1.0, 1.0); angle = std::acos(dP); if((angle * radToDeg) < angleTolerance) { @@ -110,10 +92,11 @@ bool check_for_burgers(const QuatD& betaQuat, const QuatD& alphaQuat, float64 an { return true; } - b[0] = -0.5 * gAlphaT[0][0] - 0.866025 * gAlphaT[0][1]; - b[1] = -0.5 * gAlphaT[1][0] - 0.866025 * gAlphaT[1][1]; - b[2] = -0.5 * gAlphaT[2][0] - 0.866025 * gAlphaT[2][1]; - dP = GeometryMath::CosThetaBetweenVectors(a, b); + b[0] = -0.5 * gAlphaT[0] - 0.866025 * gAlphaT[1]; + b[1] = -0.5 * gAlphaT[3] - 0.866025 * gAlphaT[4]; + b[2] = -0.5 * gAlphaT[6] - 0.866025 * gAlphaT[7]; + dP = a.cosTheta(b); + dP = std::clamp(dP, -1.0, 1.0); angle = std::acos(dP); if((angle * radToDeg) < angleTolerance) { @@ -135,7 +118,7 @@ MergeColonies::MergeColonies(DataStructure& dataStructure, const IFilter::Messag , m_InputValues(inputValues) , m_ShouldCancel(shouldCancel) , m_MessageHandler(mesgHandler) -, m_OrientationOps(LaueOps::GetAllOrientationOps()) +, m_OrientationOps(ebsdlib::LaueOps::GetAllOrientationOps()) , m_FeatureParentIds(dataStructure.getDataRefAs(inputValues->FeatureParentIdsPath)) , m_FeaturePhases(dataStructure.getDataRefAs(inputValues->FeaturePhasesPath)) , m_AvgQuats(dataStructure.getDataRefAs(inputValues->AvgQuatsPath)) @@ -148,12 +131,6 @@ MergeColonies::MergeColonies(DataStructure& dataStructure, const IFilter::Messag // ----------------------------------------------------------------------------- MergeColonies::~MergeColonies() noexcept = default; -// ----------------------------------------------------------------------------- -const std::atomic_bool& MergeColonies::getCancel() -{ - return m_ShouldCancel; -} - // ----------------------------------------------------------------------------- bool MergeColonies::growPatch(int32_t currentPatch) { @@ -167,70 +144,85 @@ bool MergeColonies::growGrouping(int32_t referenceFeature, int32_t neighborFeatu } // ----------------------------------------------------------------------------- -void MergeColonies::execute() +Result<> MergeColonies::execute() { - NeighborList& neighborlist = m_DataStructure.getDataRefAs>(m_InputValues->ContiguousNeighborListArrayPath); - NeighborList* nonContigNeighList = nullptr; + MessageHelper messageHelper(m_MessageHandler); + ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger(); + + NeighborList& featureNeighborListRef = m_DataStructure.getDataRefAs>(m_InputValues->ContiguousNeighborListArrayPath); + NeighborList* nonContigNeighListPtr = nullptr; if(m_InputValues->UseNonContiguousNeighbors) { - nonContigNeighList = m_DataStructure.getDataAs>(m_InputValues->NonContiguousNeighborListArrayPath); + nonContigNeighListPtr = m_DataStructure.getDataAs>(m_InputValues->NonContiguousNeighborListArrayPath); + } + if(nullptr == nonContigNeighListPtr) + { + return MakeErrorResult(-99342, "There was an error getting the Non-contiguous neighborlist from the DataStructure"); } - std::vector grouplist; + std::vector groupList; - int32 parentcount = 0; - int32 seed = 0; - int32 list1size = 0, list2size = 0, listsize = 0; - int32 neigh = 0; + int32 parentCount = 0; + int32 featureSeed = 0; + int32 featureNeighborListSize = 0, nonContigNeighListSize = 0; bool patchGrouping = false; - while(seed >= 0) + while(featureSeed >= 0) { - parentcount++; - seed = getSeed(parentcount); - if(seed >= 0) + parentCount++; + featureSeed = getSeed(parentCount); + if(featureSeed >= 0) { - grouplist.push_back(seed); - for(std::vector::size_type j = 0; j < grouplist.size(); j++) + groupList.push_back(featureSeed); + // Loop over the current `groupList` vector which can grow during the looping + for(std::vector::size_type j = 0; j < groupList.size(); j++) { - int32 firstfeature = grouplist[j]; - list1size = int32(neighborlist[firstfeature].size()); + int32 firstFeature = groupList[j]; + featureNeighborListSize = static_cast(featureNeighborListRef[firstFeature].size()); if(m_InputValues->UseNonContiguousNeighbors) { - list2size = nonContigNeighList->getListSize(firstfeature); + nonContigNeighListSize = nonContigNeighListPtr->getListSize(firstFeature); } + + // There are 2 kinds of NeighborLists so we loop on both of them + // k=0: FeatureNeighborList + // k=1: FeatureNeighborHood (non-contiguous neighbors) for(int32 k = 0; k < 2; k++) { + // If we are patchGrouping, then skip the first group if(patchGrouping) { k = 1; } + int32 currentListSize = 0; if(k == 0) { - listsize = list1size; + currentListSize = featureNeighborListSize; } else if(k == 1) { - listsize = list2size; + currentListSize = nonContigNeighListSize; } - for(int32 l = 0; l < listsize; l++) + // Loop over which ever NeighborList we are currently using + for(int32 l = 0; l < currentListSize; l++) { - if(k == 0) + int32 neigh = -1; + if(k == 0) // Feature Neighbor List { - neigh = neighborlist[firstfeature][l]; + neigh = featureNeighborListRef[firstFeature][l]; } - else if(k == 1) + else if(k == 1 && m_InputValues->UseNonContiguousNeighbors) // Feature NeighborHood (non-contiguous) { bool ok = false; - neigh = nonContigNeighList->getValue(firstfeature, l, ok); + neigh = nonContigNeighListPtr->getValue(firstFeature, l, ok); } - if(neigh != firstfeature) + if(neigh >= 0 && neigh != firstFeature) { - if(determineGrouping(firstfeature, neigh, parentcount)) + if(determineGrouping(firstFeature, neigh, parentCount)) { if(!patchGrouping) { - grouplist.push_back(neigh); + groupList.push_back(neigh); } } } @@ -239,35 +231,47 @@ void MergeColonies::execute() } if(patchGrouping) { - if(growPatch(parentcount)) + if(growPatch(parentCount)) { - for(std::vector::size_type j = 0; j < grouplist.size(); j++) + for(std::vector::size_type j = 0; j < groupList.size(); j++) { - int32_t firstfeature = grouplist[j]; - listsize = int32_t(neighborlist[firstfeature].size()); - for(int32_t l = 0; l < listsize; l++) + int32_t firstfeature = groupList[j]; + int32 currentListSize = static_cast(featureNeighborListRef[firstfeature].size()); + for(int32_t l = 0; l < currentListSize; l++) { - neigh = neighborlist[firstfeature][l]; + int32 neigh = featureNeighborListRef[firstfeature][l]; if(neigh != firstfeature) { - if(growGrouping(firstfeature, neigh, parentcount)) + if(growGrouping(firstfeature, neigh, parentCount)) { - grouplist.push_back(neigh); + groupList.push_back(neigh); } } } } } } + + throttledMessenger.sendThrottledMessage([&]() { return fmt::format("Parent Count: {}", parentCount); }); } - grouplist.clear(); + groupList.clear(); } + return {}; } // ----------------------------------------------------------------------------- Result<> MergeColonies::operator()() { - execute(); + // Initialize the random number generator + m_Generator = std::mt19937_64(m_InputValues->SeedValue); + m_Distribution = std::uniform_real_distribution(0.0f, 1.0f); + + // The main algorithm is in the 'execute()' method + Result<> result = execute(); + if(result.invalid()) + { + return result; + } auto active = m_DataStructure.getDataRefAs(m_InputValues->ActivePath); @@ -335,7 +339,7 @@ Result<> MergeColonies::operator()() } m_MessageHandler({IFilter::Message::Type::Info, "Adjusting Feature Ids Array...."}); - // Now adjust all the Feature Id values for each Voxel + // Now adjust all the FeatureId values for each Voxel for(usize i = 0; i < totalPoints; ++i) { cellParentIds[i] = pid[cellParentIds[i]]; @@ -351,114 +355,127 @@ Result<> MergeColonies::operator()() // ----------------------------------------------------------------------------- int32 MergeColonies::getSeed(int32 newFid) { - usize numFeatures = m_FeaturePhases.getNumberOfTuples(); + const usize numFeatures = m_FeaturePhases.getNumberOfTuples(); - std::mt19937 generator(m_InputValues->SeedValue); // Standard mersenne_twister_engine seeded - std::uniform_real_distribution distribution(0, 1); - int32 seed = -1; - int32 randFeature = 0; + int32 featureIdSeed = -1; // Precalculate some constants - usize totalFMinus1 = numFeatures - 1; + const int32 totalFMinus1 = static_cast(numFeatures) - 1; usize counter = 0; - randFeature = int32(distribution(generator) * float32(totalFMinus1)); - while(seed == -1 && counter < numFeatures) + // This section finds a feature id that has not been grouped yet. It starts by + // randomly selecting a feature id between 0 and numFeatures-1. We then start + // looping. If the initial random value is valid then we exit the loop after + // a single iteration. If that feature has already been grouped, then we add one + // to the `randFeature` value and try again. If we get to the end of the range of + // featureIds then the algorithm will loop back to featureId = 0 and start incrementing + // from there. This is reasonably efficient as we only generate random numbers + // as needed. + auto randFeature = static_cast(m_Distribution(m_Generator) * static_cast(totalFMinus1)); + while(featureIdSeed == -1 && counter < numFeatures) { if(randFeature > totalFMinus1) { randFeature = randFeature - numFeatures; } - if(m_FeatureParentIds[randFeature] == -1) + if(m_FeatureParentIds.getValue(randFeature) == -1) { - seed = randFeature; + featureIdSeed = randFeature; } randFeature++; counter++; } - if(seed >= 0) + + // Resize the created Feature Attribute Matrix + if(featureIdSeed >= 0) { - m_FeatureParentIds[seed] = newFid; - std::vector tDims(1, newFid + 1); + m_FeatureParentIds[featureIdSeed] = newFid; + const std::vector tDims = {static_cast(newFid + 1)}; m_DataStructure.getDataRefAs(m_InputValues->CellFeatureAMPath).resizeTuples(tDims); } - return seed; + return featureIdSeed; } // ----------------------------------------------------------------------------- -bool MergeColonies::determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid) +bool MergeColonies::determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid) const { - float64 w = std::numeric_limits::max(); - bool colony = false; - + // The phase is valid for both features and the neighbor feature has not been grouped yet. if(m_FeatureParentIds[neighborFeature] == -1 && m_FeaturePhases[referenceFeature] > 0 && m_FeaturePhases[neighborFeature] > 0) { usize avgQuatIdx = referenceFeature * 4; - QuatD q1(m_AvgQuats[avgQuatIdx], m_AvgQuats[avgQuatIdx + 1], m_AvgQuats[avgQuatIdx + 2], m_AvgQuats[avgQuatIdx + 3]); + const ebsdlib::QuatD q1(m_AvgQuats[avgQuatIdx], m_AvgQuats[avgQuatIdx + 1], m_AvgQuats[avgQuatIdx + 2], m_AvgQuats[avgQuatIdx + 3]); avgQuatIdx = neighborFeature * 4; - QuatD q2(m_AvgQuats[avgQuatIdx], m_AvgQuats[avgQuatIdx + 1], m_AvgQuats[avgQuatIdx + 2], m_AvgQuats[avgQuatIdx + 3]); + const ebsdlib::QuatD q2(m_AvgQuats[avgQuatIdx], m_AvgQuats[avgQuatIdx + 1], m_AvgQuats[avgQuatIdx + 2], m_AvgQuats[avgQuatIdx + 3]); - uint32 laueClass1 = m_CrystalStructures[m_FeaturePhases[referenceFeature]]; - uint32 laueClass2 = m_CrystalStructures[m_FeaturePhases[neighborFeature]]; - if(laueClass1 == laueClass2 && (laueClass1 == EbsdLib::CrystalStructure::Hexagonal_High)) + // Make sure both features are of the same Laue class and the Laue class is hexagonal + const uint32 laueClass1 = m_CrystalStructures[m_FeaturePhases[referenceFeature]]; + const uint32 laueClass2 = m_CrystalStructures[m_FeaturePhases[neighborFeature]]; + if(laueClass1 == laueClass2 && (laueClass1 == ebsdlib::CrystalStructure::Hexagonal_High)) { - OrientationD ax = m_OrientationOps[laueClass1]->calculateMisorientation(q1, q2); + ebsdlib::AxisAngleDType ax = m_OrientationOps[laueClass1]->calculateMisorientation(q1, q2); - auto rod = OrientationTransformation::ax2ro(ax); + ebsdlib::Rodrigues rod = ax.toRodrigues(); rod = m_OrientationOps[laueClass1]->getMDFFZRod(rod); - ax = OrientationTransformation::ro2ax(rod); + ax = rod.toAxisAngle(); + const float32 w = ax[3] * (Constants::k_180OverPiD); // Convert to degrees - w = ax[3] * (Constants::k_180OverPiD); - float angdiff1 = std::fabs(w - 10.53f); - float axisdiff1 = std::acos( - /*std::fabs(n1) * 0.0000f + std::fabs(n2) * 0.0000f +*/ std::fabs(ax[2]) /* * 1.0000f */); + // Test each of the special Axis-Angle relationships + // c = <0001> + float angdiff1 = std::fabs(w - 10.529f); + float axisdiff1 = std::acos(std::fabs(ax[2])); if(angdiff1 < m_AngleTolerance && axisdiff1 < m_AxisToleranceRad) { - colony = true; - } - float angdiff2 = std::fabs(w - 90.00f); - float axisdiff2 = std::acos(std::fabs(ax[0]) * 0.9958f + std::fabs(ax[1]) * 0.0917f /* + std::fabs(n3) * 0.0000f */); - if(angdiff2 < m_AngleTolerance && axisdiff2 < m_AxisToleranceRad) - { - colony = true; + m_FeatureParentIds[neighborFeature] = newFid; + return true; } + + // a2 = <-12-10> float angdiff3 = std::fabs(w - 60.00f); - float axisdiff3 = std::acos(std::fabs(ax[0]) /* * 1.0000f + std::fabs(n2) * 0.0000f + std::fabs(n3) * 0.0000f*/); + float axisdiff3 = std::acos(std::fabs(ax[0])); if(angdiff3 < m_AngleTolerance && axisdiff3 < m_AxisToleranceRad) { - colony = true; + m_FeatureParentIds[neighborFeature] = newFid; + return true; } + + // d1 at 80.97 degrees from c in the plane of (d3,c) float angdiff4 = std::fabs(w - 60.83f); float axisdiff4 = std::acos(std::fabs(ax[0]) * 0.9834f + std::fabs(ax[1]) * 0.0905f + std::fabs(ax[2]) * 0.1570f); if(angdiff4 < m_AngleTolerance && axisdiff4 < m_AxisToleranceRad) { - colony = true; + m_FeatureParentIds[neighborFeature] = newFid; + return true; } + + // d2 at 72.73 degrees from c in the plane of (a2,c) float angdiff5 = std::fabs(w - 63.26f); - float axisdiff5 = std::acos(std::fabs(ax[0]) * 0.9549f /* + std::fabs(n2) * 0.0000f */ + std::fabs(ax[2]) * 0.2969f); + float axisdiff5 = std::acos(std::fabs(ax[0]) * 0.9549f + std::fabs(ax[2]) * 0.2969f); if(angdiff5 < m_AngleTolerance && axisdiff5 < m_AxisToleranceRad) { - colony = true; + m_FeatureParentIds[neighborFeature] = newFid; + return true; } - if(colony) + + // d3 at 5.26 degrees from a2 in the basal plane + float angdiff2 = std::fabs(w - 90.00f); + float axisdiff2 = std::acos(std::fabs(ax[0]) * 0.9958f + std::fabs(ax[1]) * 0.0917f); + if(angdiff2 < m_AngleTolerance && axisdiff2 < m_AxisToleranceRad) { m_FeatureParentIds[neighborFeature] = newFid; return true; } } - else if(EbsdLib::CrystalStructure::Cubic_High == laueClass2 && EbsdLib::CrystalStructure::Hexagonal_High == laueClass1) + else if(ebsdlib::CrystalStructure::Cubic_High == laueClass2 && ebsdlib::CrystalStructure::Hexagonal_High == laueClass1) { - colony = check_for_burgers(q2, q1, m_AngleTolerance); - if(colony) + if(check_for_burgers(q2, q1, m_AngleTolerance)) { m_FeatureParentIds[neighborFeature] = newFid; return true; } } - else if(EbsdLib::CrystalStructure::Cubic_High == laueClass1 && EbsdLib::CrystalStructure::Hexagonal_High == laueClass2) + else if(ebsdlib::CrystalStructure::Cubic_High == laueClass1 && ebsdlib::CrystalStructure::Hexagonal_High == laueClass2) { - colony = check_for_burgers(q1, q2, m_AngleTolerance); - if(colony) + if(check_for_burgers(q1, q2, m_AngleTolerance)) { m_FeatureParentIds[neighborFeature] = newFid; return true; diff --git a/src/SimplnxReview/Filters/Algorithms/MergeColonies.hpp b/src/SimplnxReview/Filters/Algorithms/MergeColonies.hpp index 50563e0..6ba4b7f 100644 --- a/src/SimplnxReview/Filters/Algorithms/MergeColonies.hpp +++ b/src/SimplnxReview/Filters/Algorithms/MergeColonies.hpp @@ -13,6 +13,8 @@ #include "EbsdLib/LaueOps/LaueOps.h" +#include + namespace nx::core { struct SIMPLNXREVIEW_EXPORT MergeColoniesInputValues @@ -43,7 +45,7 @@ struct SIMPLNXREVIEW_EXPORT MergeColoniesInputValues class SIMPLNXREVIEW_EXPORT MergeColonies { - using LaueOpsShPtrType = std::shared_ptr; + using LaueOpsShPtrType = std::shared_ptr; using LaueOpsContainer = std::vector; public: @@ -57,12 +59,10 @@ class SIMPLNXREVIEW_EXPORT MergeColonies Result<> operator()(); - const std::atomic_bool& getCancel(); - protected: int getSeed(int32 newFid); - bool determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid); - void execute(); + bool determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid) const; + Result<> execute(); bool growPatch(int32 currentPatch); bool growGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid); void characterize_colonies(); @@ -80,5 +80,8 @@ class SIMPLNXREVIEW_EXPORT MergeColonies UInt32Array& m_CrystalStructures; float32 m_AxisToleranceRad = 0.0; float32 m_AngleTolerance = 1.0; + + std::mt19937_64 m_Generator = {}; + std::uniform_real_distribution m_Distribution = {}; }; } // namespace nx::core diff --git a/src/SimplnxReview/Filters/ComputeGroupingDensityFilter.cpp b/src/SimplnxReview/Filters/ComputeGroupingDensityFilter.cpp index cb5dc3a..a2f79ff 100644 --- a/src/SimplnxReview/Filters/ComputeGroupingDensityFilter.cpp +++ b/src/SimplnxReview/Filters/ComputeGroupingDensityFilter.cpp @@ -61,12 +61,17 @@ Parameters ComputeGroupingDensityFilter::parameters() const params.insertSeparator(Parameters::Separator{"Input Parameter(s)"}); params.insertLinkableParameter(std::make_unique(k_FindCheckedFeatures_Key, "Find Checked Features", "Find checked features", false)); + params.insertSeparator(Parameters::Separator{"Non-Contiguous Neighborhood Option"}); + params.insertLinkableParameter(std::make_unique(k_UseNonContiguousNeighbors_Key, "Use Non-Contiguous Neighbors", "Use non-contiguous neighborhoods for computations", false)); + params.insert(std::make_unique(k_NonContiguousNeighborListArrayPath_Key, "Non-Contiguous Neighborhoods", "List of non-contiguous neighbors for each Feature.", + DataPath{}, NeighborListSelectionParameter::AllowedTypes{DataType::int32})); + params.insertSeparator(Parameters::Separator{"Input Cell Data"}); params.insert(std::make_unique(k_ParentIdsPath_Key, "Parent Ids", "Input Cell level ParentIds", DataPath{}, ArraySelectionParameter::AllowedTypes{DataType::int32}, ArraySelectionParameter::AllowedComponentShapes{{1}})); params.insertSeparator(Parameters::Separator{"Input Feature Data"}); - params.insert(std::make_unique(k_VolumesArrayPath_Key, "Volumes", "Input Feature Volumes Data Array", DataPath{}, + params.insert(std::make_unique(k_VolumesArrayPath_Key, "Volumes", "The Feature Volumes Data Array", DataPath{}, ArraySelectionParameter::AllowedTypes{nx::core::DataType::float32}, ArraySelectionParameter::AllowedComponentShapes{{1}})); params.insert(std::make_unique(k_ContiguousNeighborListArrayPath_Key, "Contiguous Neighbor List", "List of contiguous neighbors for each Feature.", DataPath{}, @@ -75,11 +80,6 @@ Parameters ComputeGroupingDensityFilter::parameters() const params.insert(std::make_unique(k_ParentVolumesPath_Key, "Parent Volumes", "Input feature level parent volume data array", DataPath{}, ArraySelectionParameter::AllowedTypes{DataType::float32}, ArraySelectionParameter::AllowedComponentShapes{{1}})); - params.insertSeparator(Parameters::Separator{"Non-Contiguous Neighborhood Option"}); - params.insertLinkableParameter(std::make_unique(k_UseNonContiguousNeighbors_Key, "Use Non-Contiguous Neighbors", "Use non-contiguous neighborhoods for computations", false)); - params.insert(std::make_unique(k_NonContiguousNeighborListArrayPath_Key, "Non-Contiguous Neighborhoods", "Input feature level Non-Contiguous neighborhoods", - DataPath{}, NeighborListSelectionParameter::AllowedTypes{DataType::int32})); - params.insertSeparator(Parameters::Separator{"Output Feature Data"}); params.insert(std::make_unique(k_CheckedFeaturesName_Key, "Checked Features Name", "Output feature level data array to hold 'Checked Features' values", "Checked Features")); @@ -251,5 +251,4 @@ Result ComputeGroupingDensityFilter::FromSIMPLJson(const nlohmann::js return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/ComputeLocalAverageCAxisMisalignmentsFilter.cpp b/src/SimplnxReview/Filters/ComputeLocalAverageCAxisMisalignmentsFilter.cpp index a78332a..1a71d7b 100644 --- a/src/SimplnxReview/Filters/ComputeLocalAverageCAxisMisalignmentsFilter.cpp +++ b/src/SimplnxReview/Filters/ComputeLocalAverageCAxisMisalignmentsFilter.cpp @@ -295,5 +295,4 @@ Result ComputeLocalAverageCAxisMisalignmentsFilter::FromSIMPLJson(con return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/ComputeMicroTextureRegionsFilter.cpp b/src/SimplnxReview/Filters/ComputeMicroTextureRegionsFilter.cpp index e6d1993..49a4de5 100644 --- a/src/SimplnxReview/Filters/ComputeMicroTextureRegionsFilter.cpp +++ b/src/SimplnxReview/Filters/ComputeMicroTextureRegionsFilter.cpp @@ -159,5 +159,4 @@ Result ComputeMicroTextureRegionsFilter::FromSIMPLJson(const nlohmann return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/ComputeSaltykovSizesFilter.cpp b/src/SimplnxReview/Filters/ComputeSaltykovSizesFilter.cpp index 868ebbd..2f7372f 100644 --- a/src/SimplnxReview/Filters/ComputeSaltykovSizesFilter.cpp +++ b/src/SimplnxReview/Filters/ComputeSaltykovSizesFilter.cpp @@ -147,5 +147,4 @@ Result ComputeSaltykovSizesFilter::FromSIMPLJson(const nlohmann::json return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/GroupMicroTextureRegionsFilter.cpp b/src/SimplnxReview/Filters/GroupMicroTextureRegionsFilter.cpp index 0266a03..d8cb659 100644 --- a/src/SimplnxReview/Filters/GroupMicroTextureRegionsFilter.cpp +++ b/src/SimplnxReview/Filters/GroupMicroTextureRegionsFilter.cpp @@ -56,17 +56,21 @@ Parameters GroupMicroTextureRegionsFilter::parameters() const // Create the parameter descriptors that are needed for this filter params.insertSeparator(Parameters::Separator{"Input Parameter(s)"}); + + params.insert(std::make_unique(k_UseRunningAverage_Key, "Group C-Axes With Running Average", "Group C-Axes With Running Average", true)); + params.insert(std::make_unique(k_CAxisTolerance_Key, "C-Axis Alignment Tolerance (Degrees)", "C-Axis Alignment Tolerance (Degrees)", 0.0f)); + params.insert(std::make_unique(k_ContiguousNeighborListArrayPath_Key, "Contiguous Neighbor List", "List of contiguous neighbors for each Feature.", DataPath{}, + NeighborListSelectionParameter::AllowedTypes{DataType::int32})); + + params.insertSeparator(Parameters::Separator{"Non-Contiguous Neighborhood Option"}); params.insertLinkableParameter(std::make_unique(k_UseNonContiguousNeighbors_Key, "Use Non-Contiguous Neighbors", "Use non-contiguous neighborhoods", false)); params.insert(std::make_unique(k_NonContiguousNeighborListArrayPath_Key, "Non-Contiguous Neighbor List", "List of non-contiguous neighbors for each Feature.", DataPath{}, NeighborListSelectionParameter::AllowedTypes{DataType::int32})); - params.insert(std::make_unique(k_ContiguousNeighborListArrayPath_Key, "Contiguous Neighbor List", "List of contiguous neighbors for each Feature.", DataPath{}, - NeighborListSelectionParameter::AllowedTypes{DataType::int32})); - params.insert(std::make_unique(k_UseRunningAverage_Key, "Group C-Axes With Running Average", "Group C-Axes With Running Average", true)); - params.insert(std::make_unique(k_CAxisTolerance_Key, "C-Axis Alignment Tolerance (Degrees)", "C-Axis Alignment Tolerance (Degrees)", 0.0f)); params.insertSeparator(Parameters::Separator{"Random Number Seed Parameters"}); params.insertLinkableParameter(std::make_unique(k_UseSeed_Key, "Use Seed for Random Generation", "When true the user will be able to put in a seed for random generation", false)); params.insert(std::make_unique>(k_SeedValue_Key, "Seed", "The seed fed into the random generator", std::mt19937::default_seed)); + params.insert(std::make_unique(k_SeedArrayName_Key, "Stored Seed Value Array Name", "Name of array holding the seed value", "_Group_MicroTexture_Regions_Seed_Value_")); params.insertSeparator(Parameters::Separator{"Input Cell Data"}); params.insert(std::make_unique(k_FeatureIdsArrayPath_Key, "Cell Feature Ids", "Data Array that specifies to which Feature each Element belongs", DataPath{}, @@ -90,8 +94,6 @@ Parameters GroupMicroTextureRegionsFilter::parameters() const params.insert(std::make_unique(k_NewCellFeatureAttributeMatrixName_Key, "Created Microtexture Feature Attribute Matrix", "Output Feature Attribute Matrix for Microtexture Regions", DataPath{})); params.insert(std::make_unique(k_ActiveArrayName_Key, "Active Array Name", "Output Active Array", "Active")); - params.insert(std::make_unique(k_SeedArrayName_Key, "Stored Seed Value Array Name", "Output data array to store the random number seed value", - "_Group_MicroTexture_Regions_Seed_Value_")); // Associate the Linkable Parameter(s) to the children parameters that they control params.linkParameters(k_UseNonContiguousNeighbors_Key, k_NonContiguousNeighborListArrayPath_Key, true); @@ -243,5 +245,4 @@ Result GroupMicroTextureRegionsFilter::FromSIMPLJson(const nlohmann:: return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/InterpolateValuesToUnstructuredGridFilter.cpp b/src/SimplnxReview/Filters/InterpolateValuesToUnstructuredGridFilter.cpp index ae91e87..5b34607 100644 --- a/src/SimplnxReview/Filters/InterpolateValuesToUnstructuredGridFilter.cpp +++ b/src/SimplnxReview/Filters/InterpolateValuesToUnstructuredGridFilter.cpp @@ -172,5 +172,4 @@ Result InterpolateValuesToUnstructuredGridFilter::FromSIMPLJson(const return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core diff --git a/src/SimplnxReview/Filters/MergeColoniesFilter.cpp b/src/SimplnxReview/Filters/MergeColoniesFilter.cpp index 559e08c..6e00796 100644 --- a/src/SimplnxReview/Filters/MergeColoniesFilter.cpp +++ b/src/SimplnxReview/Filters/MergeColoniesFilter.cpp @@ -250,5 +250,4 @@ Result MergeColoniesFilter::FromSIMPLJson(const nlohmann::json& json) return ConvertResultTo(std::move(conversionResult), std::move(args)); } - } // namespace nx::core