diff --git a/src/Plugins/OrientationAnalysis/docs/ComputeAvgOrientationsFilter.md b/src/Plugins/OrientationAnalysis/docs/ComputeAvgOrientationsFilter.md index f2870baaaf..b3e261622d 100644 --- a/src/Plugins/OrientationAnalysis/docs/ComputeAvgOrientationsFilter.md +++ b/src/Plugins/OrientationAnalysis/docs/ComputeAvgOrientationsFilter.md @@ -6,16 +6,67 @@ Statistics (Crystallography) ## Description -This **Filter** determines the average orientation of each **Feature** by the following algorithm: +This **Filter** computes the average orientation of each **Feature** using one or more of three available averaging methods. Each method can be independently enabled, and their results are stored in separate output arrays. -1. Gather all **Elements** that belong to the **Feature** -2. Using the symmetry operators of the phase of the **Feature**, rotate the quaternion of the **Feature**'s first **Element** into the *Fundamental Zone* nearest to the origin -3. Rotate each subsequent **Elements**'s quaternion (with same symmetry operators) looking for the quaternion closest to the quaternion selected in Step 2 -4. Average the rotated quaternions for all **Elements** and store as the average for the **Feature** +### Method 1: Rodrigues Average (Original) -*Note:* The process of finding the nearest quaternion in Step 3 is to account for the periodicity of orientation space, which would cause problems in the averaging if all quaternions were forced to be rotated into the same *Fundamental Zone* +This is the original averaging algorithm. It determines the average orientation of each **Feature** by: -*Note:* The quaternions can be averaged with a simple average because the quaternion space is not distorted like Euler space. +1. Gathering all **Elements** that belong to the **Feature** +2. Using the symmetry operators of the phase of the **Feature**, rotating the quaternion of the **Feature**'s first **Element** into the *Fundamental Zone* nearest to the origin +3. Rotating each subsequent **Element**'s quaternion (with same symmetry operators) looking for the quaternion closest to the current running average +4. Accumulating a running sum of the nearest quaternions +5. Dividing the accumulated quaternion sum by the count and normalizing to produce the average + +The process of finding the nearest quaternion in Step 3 accounts for the periodicity of orientation space, which would cause problems in the averaging if all quaternions were forced to be rotated into the same *Fundamental Zone*. The quaternions can be averaged with a simple summation-based average because quaternion space is not distorted like Euler space. + +**Outputs:** Average Quaternions, Average Euler Angles (Bunge convention Z-X-Z) + +### Method 2: Von Mises-Fisher (vMF) Average + +The von Mises-Fisher distribution is a probability distribution on the surface of a unit hypersphere in *p*-dimensional space. For orientation averaging, the relevant case is the unit quaternion sphere (*p* = 4). The vMF distribution is parameterized by: + +- **mu** (mean direction): A unit quaternion representing the central tendency of the distribution. This is the estimated average orientation. +- **kappa** (concentration parameter): A non-negative scalar that characterizes how tightly the orientations are clustered around the mean. Intuitively, kappa is to the vMF distribution what the full-width-at-half-maximum (FWHM) is to a Gaussian distribution: it is a measure of how narrow or tight the distribution is. Higher kappa values indicate tighter clustering (less spread); kappa = 0 corresponds to a uniform distribution on the sphere. + +The vMF probability density for a unit vector **x** given mean direction **mu** and concentration **kappa** is proportional to exp(kappa * mu^T * x). This makes it the spherical analogue of the Gaussian distribution on a flat space. + +The filter estimates the vMF parameters using an **Expectation-Maximization (EM)** algorithm. All element quaternions belonging to a feature are first reduced to the *Fundamental Zone* using the crystal symmetry operators. The EM procedure then iteratively refines the estimates of **mu** (the average orientation quaternion) and **kappa** (the concentration). + +**Outputs:** Average Quaternions, Average Euler Angles (Bunge convention Z-X-Z), Kappa Values + +### Method 3: Watson Average + +The Watson distribution is a probability distribution on the unit sphere that is **antipodally symmetric**, meaning it treats **x** and **-x** as equivalent. This property makes it particularly well-suited for orientation data represented as quaternions, since quaternions **q** and **-q** represent the same physical rotation. + +The Watson distribution is parameterized by: + +- **mu** (mean axis): A unit quaternion representing the principal axis of the distribution. This is the estimated average orientation. +- **kappa** (concentration parameter): A scalar that controls the concentration of the distribution around the mean axis. As with the vMF distribution, kappa is analogous to the full-width-at-half-maximum (FWHM) of a Gaussian: it measures how narrow or tight the distribution is. For positive kappa, the distribution is bipolar (concentrated around +/-mu); for negative kappa, it is girdle-shaped (concentrated in the great circle perpendicular to mu). + +The Watson probability density for a unit vector **x** is proportional to exp(kappa * (mu^T * x)^2). The key difference from the von Mises-Fisher distribution is the squared dot product, which enforces the antipodal symmetry. + +Like the vMF method, the filter estimates Watson parameters using an **Expectation-Maximization (EM)** algorithm operating on fundamental-zone-reduced quaternions. + +**Outputs:** Average Quaternions, Average Euler Angles (Bunge convention Z-X-Z), Kappa Values + +### Hard-Coded Algorithm Parameters + +The following parameters are currently hard-coded in the implementation and are not user-configurable: + +| Parameter | Value | Description | +|-----------|-------|-------------| +| **Random Seed** | 43514 | Seed for the random number generator used in the EM algorithm. Because this is fixed, the vMF and Watson results are deterministic across runs. | +| **EM Iterations** | 5 | Number of Expectation-Maximization outer iterations. Controls how many times the full EM cycle is repeated. | +| **Iterations** | 10 | Number of inner iterations per EM cycle. Controls the refinement within each EM step. | + +These values may be exposed as user-configurable parameters in a future release. + +### Special Cases + +- **Features with a single element:** For the vMF and Watson methods, if a feature contains only one element orientation, the EM algorithm is skipped entirely and the single quaternion is used directly as the average. The kappa value is set to 0 in this case. +- **Features with zero elements:** Features with no elements (phase <= 0 for all voxels) will have their output arrays initialized to NaN (for vMF/Watson) or identity quaternion / zero Euler angles (for Rodrigues). +- **Phase indexing:** The filter requires that phase values be > 0 for elements to be included in the averaging. Phase index 0 is reserved for "Unknown" in the Crystal Structures array and is always skipped. % Auto generated parameter table will be inserted here diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.cpp index fcf27bb5bd..40b9b53e0c 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.cpp @@ -2,7 +2,10 @@ #include "simplnx/DataStructure/AttributeMatrix.hpp" #include "simplnx/Utilities/DataArrayUtilities.hpp" +#include "simplnx/Utilities/MessageHelper.hpp" +#include "simplnx/Utilities/ParallelDataAlgorithm.hpp" +#include #include #include @@ -12,17 +15,167 @@ using namespace nx::core; namespace { -// std::ostream& operator<<(std::ostream& os, const ebsdlib::QuatF& q) -// { -// os << ", " << q.x() << ", " << q.y() << ", " << q.z() << ", " << q.w(); -// return os; -// } -// -// std::ostream& operator<<(std::ostream& os, const ebsdlib::QuatD& q) -// { -// os << ", " << q.x() << ", " << q.y() << ", " << q.z() << ", " << q.w(); -// return os; -// } +class VmfWatsonSamplingImpl +{ +public: + VmfWatsonSamplingImpl(ComputeAvgOrientations* filter, const ComputeAvgOrientationsInputValues* inputPtr, DataStructure& dataStruture, const std::vector& featureNumVoxels, + const std::map& featureIdToPhaseMap) + : m_Filter(filter) + , m_InputValues(inputPtr) + , m_DataStructure(dataStruture) + , m_FeatureNumVoxels(featureNumVoxels) + , m_FeatureIdToPhaseMap(featureIdToPhaseMap) + { + } + + virtual ~VmfWatsonSamplingImpl() = default; + + void operator()(const Range& range) const + { + // Input FeatureIds + Input Orientations. All these should come from the same Attribute Matrix or have the same number of tuples + Int32AbstractDataStore& featureIdsRef = m_DataStructure.getDataRefAs(m_InputValues->cellFeatureIdsArrayPath).getDataStoreRef(); + // Int32AbstractDataStore& phasesRef = m_DataStructure.getDataRefAs(m_InputValues->cellPhasesArrayPath).getDataStoreRef(); + Float32AbstractDataStore& quatsRef = m_DataStructure.getDataRefAs(m_InputValues->cellQuatsArrayPath).getDataStoreRef(); + // Ensemble Level Data + UInt32AbstractDataStore& xtalRef = m_DataStructure.getDataRefAs(m_InputValues->crystalStructuresArrayPath).getDataStoreRef(); + + // Output vMF Data + Float32AbstractDataStore* vmfQuatPtr = nullptr; + Float32AbstractDataStore* vmfEulerPtr = nullptr; + Float32AbstractDataStore* vmfKappaPtr = nullptr; + if(m_InputValues->useVonMisesAverage) + { + vmfQuatPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFQuatsArrayPath).getDataStorePtr().lock().get(); + vmfEulerPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFEulerAnglesArrayPath).getDataStorePtr().lock().get(); + vmfKappaPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFKappaArrayPath).getDataStorePtr().lock().get(); + } + + // Output Watson Data + Float32AbstractDataStore* watsonQuatPtr = nullptr; + Float32AbstractDataStore* watsonEulerPtr = nullptr; + Float32AbstractDataStore* watsonKappaPtr = nullptr; + if(m_InputValues->useWatsonAverage) + { + watsonQuatPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonQuatsArrayPath).getDataStorePtr().lock().get(); + watsonEulerPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonEulerAnglesArrayPath).getDataStorePtr().lock().get(); + watsonKappaPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonKappaArrayPath).getDataStorePtr().lock().get(); + } + + usize numVoxels = featureIdsRef.getNumberOfTuples(); + + std::vector ops = ebsdlib::LaueOps::GetAllOrientationOps(); + + std::vector fzQuats; + for(usize featureId = range.min(); featureId < range.max(); featureId++) + { + // If the size is 0 then skip to the next feature + if(m_FeatureNumVoxels[featureId] == 0) + { + continue; + } + + const int32 phaseIdx = m_FeatureIdToPhaseMap.at(static_cast(featureId)); + const uint32 laueClass = xtalRef[phaseIdx]; + ebsdlib::LaueOps::Pointer op = ops[laueClass]; + + fzQuats.clear(); + fzQuats.reserve(m_FeatureNumVoxels[featureId]); + + // Loop over every "voxel" (although the user could just be passing in an array + // they want to find the average orientation of + for(usize voxelIdx = 0; voxelIdx < numVoxels; voxelIdx++) + { + // If the feature Id of the voxel matches the current feature Id, then grab that orientation + if(featureIdsRef[voxelIdx] == featureId) + { + const ebsdlib::QuatD q1(quatsRef[voxelIdx * 4], quatsRef[voxelIdx * 4 + 1], quatsRef[voxelIdx * 4 + 2], quatsRef[voxelIdx * 4 + 3]); + fzQuats.push_back(op->getFZQuat(q1)); // Fundamental Zone Reduction + } + } + + // Now that we have all the orientations for the given featureId we can compute the averages + if(m_InputValues->useVonMisesAverage) + { + uint32_t seed = m_InputValues->RandomSeed; // This should be a user facing options + ebsdlib::QuatD muhat = ebsdlib::QuatD::identity(); + double kappahat = 0.0; + + if(fzQuats.size() == 1) + { + muhat = fzQuats[0]; + } + else if(!fzQuats.empty()) + { + ebsdlib::DirectionalStats directionalStats("VMF", op); + int numEmIterations = m_InputValues->NumEMIterations; // At some point this should be a user-defined input + int numIterations = m_InputValues->NumIterations; // At some point this should be a user-defined input + directionalStats.setNumEM(numEmIterations); + directionalStats.setNumIter(numIterations); + directionalStats.setQuatArray(fzQuats); + directionalStats.EMforDS(seed, muhat, kappahat, false); + muhat.positiveOrientation(); + } + + vmfQuatPtr->setValue(featureId * 4, muhat.x()); + vmfQuatPtr->setValue(featureId * 4 + 1, muhat.y()); + vmfQuatPtr->setValue(featureId * 4 + 2, muhat.z()); + vmfQuatPtr->setValue(featureId * 4 + 3, muhat.w()); + + ebsdlib::EulerDType euler = muhat.toEuler(); + vmfEulerPtr->setValue(featureId * 3, euler[0]); + vmfEulerPtr->setValue(featureId * 3 + 1, euler[1]); + vmfEulerPtr->setValue(featureId * 3 + 2, euler[2]); + + vmfKappaPtr->setValue(featureId, kappahat); + } + + if(m_InputValues->useWatsonAverage) + { + uint32_t seed = m_InputValues->RandomSeed; // This should be a user facing options + ebsdlib::QuatD muhat = ebsdlib::QuatD::identity(); + double kappahat = 0.0; + + // Check if there is only a single orientation... + if(fzQuats.size() == 1) + { + muhat = fzQuats[0]; + } + else if(!fzQuats.empty()) + { + ebsdlib::DirectionalStats directionalStats("WAT", op); + int numEmIterations = m_InputValues->NumEMIterations; // At some point this should be a user-defined input + int numIterations = m_InputValues->NumIterations; // At some point this should be a user-defined input + directionalStats.setNumEM(numEmIterations); + directionalStats.setNumIter(numIterations); + directionalStats.setQuatArray(fzQuats); + directionalStats.EMforDS(seed, muhat, kappahat, false); + muhat.positiveOrientation(); + } + + watsonQuatPtr->setValue(featureId * 4, muhat.x()); + watsonQuatPtr->setValue(featureId * 4 + 1, muhat.y()); + watsonQuatPtr->setValue(featureId * 4 + 2, muhat.z()); + watsonQuatPtr->setValue(featureId * 4 + 3, muhat.w()); + + ebsdlib::EulerDType euler = muhat.toEuler(); + watsonEulerPtr->setValue(featureId * 3, euler[0]); + watsonEulerPtr->setValue(featureId * 3 + 1, euler[1]); + watsonEulerPtr->setValue(featureId * 3 + 2, euler[2]); + + watsonKappaPtr->setValue(featureId, kappahat); + } + + m_Filter->sendThreadSafeProgressMessage(1); + } + } + +private: + ComputeAvgOrientations* m_Filter = nullptr; + const ComputeAvgOrientationsInputValues* m_InputValues = nullptr; + DataStructure& m_DataStructure; + const std::vector& m_FeatureNumVoxels; + const std::map& m_FeatureIdToPhaseMap; +}; template void UpdateQuaternionArray(AbstractDataStore& quatArray, const ebsdlib::Quaternion& quat, int32 tupleIndex) @@ -56,8 +209,162 @@ ComputeAvgOrientations::ComputeAvgOrientations(DataStructure& dataStructure, con // ----------------------------------------------------------------------------- ComputeAvgOrientations::~ComputeAvgOrientations() noexcept = default; +// ----------------------------------------------------------------------------- +void ComputeAvgOrientations::sendThreadSafeProgressMessage(usize counter) +{ + std::lock_guard guard(m_ProgressMessage_Mutex); + + m_ProgressCounter += counter; + auto now = std::chrono::steady_clock::now(); + if(std::chrono::duration_cast(now - m_InitialPoint).count() < 1000) + { + return; + } + + auto progressInt = static_cast((static_cast(m_ProgressCounter) / static_cast(m_NumberOfFeatures)) * 100.0f); + std::string ss = fmt::format("{}% Complete", progressInt); + m_MessageHandler(IFilter::Message::Type::Info, ss); + + m_LastProgressInt = progressInt; + m_InitialPoint = std::chrono::steady_clock::now(); +} + // ----------------------------------------------------------------------------- Result<> ComputeAvgOrientations::operator()() +{ + const auto& featureIds = m_DataStructure.getDataRefAs(m_InputValues->cellFeatureIdsArrayPath); + if(m_DataStructure.containsData(m_InputValues->avgEulerAnglesArrayPath)) + { + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->avgEulerAnglesArrayPath, featureIds, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + m_NumberOfFeatures = m_DataStructure.getDataRefAs(m_InputValues->avgEulerAnglesArrayPath).getNumberOfTuples(); + } + else if(m_DataStructure.containsData(m_InputValues->VMFEulerAnglesArrayPath)) + { + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->VMFEulerAnglesArrayPath, featureIds, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + m_NumberOfFeatures = m_DataStructure.getDataRefAs(m_InputValues->VMFEulerAnglesArrayPath).getNumberOfTuples(); + } + else if(m_DataStructure.containsData(m_InputValues->WatsonEulerAnglesArrayPath)) + { + auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->WatsonEulerAnglesArrayPath, featureIds, false, m_MessageHandler); + if(validateNumFeatResult.invalid()) + { + return validateNumFeatResult; + } + m_NumberOfFeatures = m_DataStructure.getDataRefAs(m_InputValues->WatsonEulerAnglesArrayPath).getNumberOfTuples(); + } + else + { + return MakeErrorResult(-54670, "A valid Feature level array that stores results was not found."); + } + + MessageHelper messageHelper(m_MessageHandler); + + Result<> result; + if(m_InputValues->useRodriguesAverage) + { + messageHelper.sendMessage("Computing Rodrigues Average Orientations"); + + result = computeRodriguesAverage(); + if(result.invalid()) + { + return result; + } + } + if(m_InputValues->useVonMisesAverage || m_InputValues->useWatsonAverage) + { + if(m_InputValues->useVonMisesAverage && !m_InputValues->useWatsonAverage) + { + messageHelper.sendMessage("Computing von-Mises Fisher Average Orientations"); + } + if(!m_InputValues->useVonMisesAverage && m_InputValues->useWatsonAverage) + { + messageHelper.sendMessage("Computing Watson Average Orientations"); + } + if(m_InputValues->useVonMisesAverage && m_InputValues->useWatsonAverage) + { + messageHelper.sendMessage("Computing von-Mises Fisher and Watson Average Orientations"); + } + + result = computeVmfWatsonAverage(); + if(result.invalid()) + { + return result; + } + } + + return {}; +} + +// ----------------------------------------------------------------------------- +Result<> ComputeAvgOrientations::computeVmfWatsonAverage() +{ + // Input Data + auto& featureIds = m_DataStructure.getDataRefAs(m_InputValues->cellFeatureIdsArrayPath); + auto& phases = m_DataStructure.getDataRefAs(m_InputValues->cellPhasesArrayPath); + + const size_t totalVoxels = featureIds.getNumberOfTuples(); + + // Run through the "voxels" and compute the number of voxels for each feature + std::vector featureNumVoxels(m_NumberOfFeatures, 0); + std::map featureIdToPhaseMap; + for(size_t i = 0; i < totalVoxels; i++) + { + const int32_t currentFeatureId = featureIds[i]; + const int32_t currentPhase = phases[i]; + if(currentPhase > 0) + { + featureNumVoxels[currentFeatureId]++; + featureIdToPhaseMap[currentFeatureId] = currentPhase; + } + } + + // Initialize the output arrays + // Output vMF Data + Float32AbstractDataStore* vmfQuatPtr = nullptr; + Float32AbstractDataStore* vmfEulerPtr = nullptr; + Float32AbstractDataStore* vmfKappaPtr = nullptr; + if(m_InputValues->useVonMisesAverage) + { + vmfQuatPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFQuatsArrayPath).getDataStorePtr().lock().get(); + vmfEulerPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFEulerAnglesArrayPath).getDataStorePtr().lock().get(); + vmfKappaPtr = m_DataStructure.getDataRefAs(m_InputValues->VMFKappaArrayPath).getDataStorePtr().lock().get(); + vmfQuatPtr->fill(std::numeric_limits::quiet_NaN()); + vmfEulerPtr->fill(std::numeric_limits::quiet_NaN()); + vmfKappaPtr->fill(std::numeric_limits::quiet_NaN()); + } + + // Output Watson Data + Float32AbstractDataStore* watsonQuatPtr = nullptr; + Float32AbstractDataStore* watsonEulerPtr = nullptr; + Float32AbstractDataStore* watsonKappaPtr = nullptr; + if(m_InputValues->useWatsonAverage) + { + watsonQuatPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonQuatsArrayPath).getDataStorePtr().lock().get(); + watsonEulerPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonEulerAnglesArrayPath).getDataStorePtr().lock().get(); + watsonKappaPtr = m_DataStructure.getDataRefAs(m_InputValues->WatsonKappaArrayPath).getDataStorePtr().lock().get(); + watsonQuatPtr->fill(std::numeric_limits::quiet_NaN()); + watsonEulerPtr->fill(std::numeric_limits::quiet_NaN()); + watsonKappaPtr->fill(std::numeric_limits::quiet_NaN()); + } + + // Allow data-based parallelization + ParallelDataAlgorithm dataAlg; + dataAlg.setRange(0, m_NumberOfFeatures); + dataAlg.execute(VmfWatsonSamplingImpl(this, m_InputValues, m_DataStructure, featureNumVoxels, featureIdToPhaseMap)); + + return {}; +} + +// ----------------------------------------------------------------------------- +Result<> ComputeAvgOrientations::computeRodriguesAverage() { std::vector orientationOps = ebsdlib::LaueOps::GetAllOrientationOps(); @@ -72,11 +379,6 @@ Result<> ComputeAvgOrientations::operator()() const size_t totalPoints = featureIds.getNumberOfTuples(); - auto validateNumFeatResult = ValidateFeatureIdsToFeatureAttributeMatrixIndexing(m_DataStructure, m_InputValues->avgQuatsArrayPath, featureIds, false, m_MessageHandler); - if(validateNumFeatResult.invalid()) - { - return validateNumFeatResult; - } size_t totalFeatures = avgQuats.getNumberOfTuples(); std::vector counts(totalFeatures, 0.0f); @@ -96,7 +398,18 @@ Result<> ComputeAvgOrientations::operator()() } const int32_t currentFeatureId = featureIds[i]; const int32_t currentPhase = phases[i]; - if(currentFeatureId > 0 && currentPhase > 0) + // As long as we have a valid `currentPhase` value which is used as an index + // into the CrystalStructures array. We ALWAYS ignore the first value in the + // CrystalStructures array. So therefore the `currentPhase` MUST be > 0. + // We can use `currentFeatureId = 0` because if someone is just wanting to compute + // the average of a bunch of orientations they may have labeled the "FeatureIds = 0" + // for all the values. The most important value is the `currentPhase` which for + // the grand majority of historical data should be > 0. + // Now in theory someone could absolutely manually import data into an "Ensemble" + // Array and NOT have the zero index as `unknown` in which case this check will + // fail them and they will not compute anything most likely. The documentation + // for the filter should be updated to cover these use-cases. + if(currentPhase > 0) { const uint32 xtal = crystalStructures[currentPhase]; counts[currentFeatureId] += 1.0f; @@ -117,7 +430,7 @@ Result<> ComputeAvgOrientations::operator()() } } - for(size_t featureId = 1; featureId < totalFeatures; featureId++) + for(size_t featureId = 0; featureId < totalFeatures; featureId++) { if(m_ShouldCancel) { @@ -127,15 +440,15 @@ Result<> ComputeAvgOrientations::operator()() if(counts[featureId] == 0.0f) { UpdateQuaternionArray(avgQuats, identityQuat, featureId); + continue; } ebsdlib::QuatF curAvgQuat(avgQuats[featureId * 4], avgQuats[featureId * 4 + 1], avgQuats[featureId * 4 + 2], avgQuats[featureId * 4 + 3]); curAvgQuat = curAvgQuat.scalarDivide(counts[featureId]); - curAvgQuat = curAvgQuat.normalize().getPositiveOrientation(); + curAvgQuat = curAvgQuat.normalize().getPositiveOrientation(); // Be sure the Quaterion is in the Northern Hemisphere UpdateQuaternionArray(avgQuats, curAvgQuat, featureId); - // Update the value for the average Euler. Be sure to make sure the Quaterion is in the northern hemisphere - // before converting it to a Euler Angle + // Update the value for the average Euler. ebsdlib::EulerFType eu = ebsdlib::QuaternionFType(curAvgQuat).toEuler(); UpdateEulerArray(avgEuler, eu, featureId); } diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.hpp index e7b58b36d6..b78f634067 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeAvgOrientations.hpp @@ -22,8 +22,24 @@ struct ORIENTATIONANALYSIS_EXPORT ComputeAvgOrientationsInputValues DataPath cellPhasesArrayPath; DataPath cellQuatsArrayPath; DataPath crystalStructuresArrayPath; + + bool useRodriguesAverage; + bool useVonMisesAverage; + bool useWatsonAverage; DataPath avgQuatsArrayPath; DataPath avgEulerAnglesArrayPath; + + DataPath VMFQuatsArrayPath; + DataPath VMFEulerAnglesArrayPath; + DataPath VMFKappaArrayPath; + + DataPath WatsonQuatsArrayPath; + DataPath WatsonEulerAnglesArrayPath; + DataPath WatsonKappaArrayPath; + + uint32 RandomSeed = 43514; + int32 NumEMIterations = 5; + int32 NumIterations = 10; }; /** @@ -42,12 +58,24 @@ class ORIENTATIONANALYSIS_EXPORT ComputeAvgOrientations Result<> operator()(); + void sendThreadSafeProgressMessage(usize counter); + protected: private: DataStructure& m_DataStructure; const IFilter::MessageHandler& m_MessageHandler; const std::atomic_bool& m_ShouldCancel; const ComputeAvgOrientationsInputValues* m_InputValues = nullptr; + + Result<> computeRodriguesAverage(); + Result<> computeVmfWatsonAverage(); + + // Thread safe Progress Message + std::chrono::steady_clock::time_point m_InitialPoint = std::chrono::steady_clock::now(); + mutable std::mutex m_ProgressMessage_Mutex; + size_t m_NumberOfFeatures = 0; + size_t m_ProgressCounter = 0; + size_t m_LastProgressInt = 0; }; } // namespace nx::core diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeFZQuaternions.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeFZQuaternions.cpp index ed21b8aef6..799591df16 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeFZQuaternions.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeFZQuaternions.cpp @@ -153,8 +153,8 @@ Result<> ComputeFZQuaternions::operator()() } else if(maskArray->getDataType() == DataType::uint8) { - UInt32Array* goodVoxelsArray = m_DataStructure.getDataAs(m_InputValues->MaskArrayPath); - dataAlg.execute(::GenerateFZQuatsImpl(quatArray, phaseArray, xtalArray, numPhases, goodVoxelsArray, fzQuatArray, m_ShouldCancel, warningCount)); + UInt8Array* goodVoxelsArray = m_DataStructure.getDataAs(m_InputValues->MaskArrayPath); + dataAlg.execute(::GenerateFZQuatsImpl(quatArray, phaseArray, xtalArray, numPhases, goodVoxelsArray, fzQuatArray, m_ShouldCancel, warningCount)); } else if(maskArray->getDataType() == DataType::int8) { diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.cpp index bbd8cd9226..4e41f086e4 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.cpp @@ -2,16 +2,31 @@ #include "OrientationAnalysis/Filters/Algorithms/ConvertOrientations.hpp" -#include - #include "simplnx/DataStructure/Geometry/VertexGeom.hpp" #include "simplnx/DataStructure/IDataArray.hpp" #include "simplnx/Utilities/DataArrayUtilities.hpp" +#include "simplnx/Utilities/FilterUtilities.hpp" + +#include #include using namespace nx::core; +namespace +{ +struct CopyDataFunctor +{ + template + void operator()(const IDataArray* srcIArray, IDataArray* destIArray) + { + const auto& srcArray = srcIArray->template getIDataStoreRefAs>(); + auto& destArray = destIArray->template getIDataStoreRefAs>(); + destArray.copyFrom(0, srcArray, 0, srcArray.getNumberOfTuples()); + } +}; +} // namespace + // ----------------------------------------------------------------------------- ConvertOrientationsToVertexGeometry::ConvertOrientationsToVertexGeometry(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ConvertOrientationsToVertexGeometryInputValues* inputValues) @@ -98,5 +113,20 @@ Result<> ConvertOrientationsToVertexGeometry::operator()() vertices.setComponent(i, 2, static_cast(st[2])); } + // Copy over the DataArrays to the new Vertex Geometry + ShapeType verticesTupleShape = vertices.getTupleShape(); + DataPath vertexAttrMatrixPath = m_InputValues->OutputVertexGeometryPath.createChildPath(m_InputValues->OutputVertexAttrMatrixName); + for(const auto& sourceDataPath : m_InputValues->DataPathCopySources) + { + auto* sourceDataArrayPtr = m_DataStructure.getDataAs(sourceDataPath); + DataPath destinationDataPath = vertexAttrMatrixPath.createChildPath(sourceDataArrayPtr->getName()); + auto* destinationDataArrayPtr = m_DataStructure.getDataAs(destinationDataPath); + ExecuteDataFunction(CopyDataFunctor{}, sourceDataArrayPtr->getDataType(), sourceDataArrayPtr, destinationDataArrayPtr); + // auto& destinationDataArray = m_DataStructure.getDataRefAs(destinationDataPath); + // This does not resize anything (at least it had better not), but is + // a round-about way to set the Tuple Shape on the destination array + destinationDataArrayPtr->resizeTuples(verticesTupleShape); + } + return {}; } diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.hpp index 5a42f3a652..6a3896dce0 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ConvertOrientationsToVertexGeometry.hpp @@ -21,7 +21,7 @@ struct ORIENTATIONANALYSIS_EXPORT ConvertOrientationsToVertexGeometryInputValues { ebsdlib::orientations::Type InputOrientationType; DataPath InputOrientationArrayPath; - std::vector CopyVertexArrayPaths; + std::vector DataPathCopySources; bool ConvertToFundamentalZone; DataPath CellPhasesArrayPath; DataPath CrystalStructuresArrayPath; diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.cpp index d1693f8edd..2e9814171f 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.cpp @@ -6,6 +6,7 @@ #include "simplnx/Filter/Actions/CreateArrayAction.hpp" #include "simplnx/Parameters/ArraySelectionParameter.hpp" #include "simplnx/Parameters/AttributeMatrixSelectionParameter.hpp" +#include "simplnx/Parameters/BoolParameter.hpp" #include "simplnx/Utilities/SIMPLConversion.hpp" @@ -51,6 +52,11 @@ Parameters ComputeAvgOrientationsFilter::parameters() const Parameters params; // Create the parameter descriptors that are needed for this filter + params.insertSeparator(Parameters::Separator{"Input Parameter(s)"}); + params.insertLinkableParameter(std::make_unique(k_UseRodriguesAverage_Key, "Compute Rodrigues Average", "The original algorithm.", true)); + params.insertLinkableParameter(std::make_unique(k_UseVonMisesFisher_Key, "Compute von Mises-Fisher Average", "The von Mises Fisher average algorithm.", false)); + params.insertLinkableParameter(std::make_unique(k_UseWatson_Key, "Compute Watson Average", "The Watson average algorithm.", false)); + 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({"Cell Data", "FeatureIds"}), ArraySelectionParameter::AllowedTypes{DataType::int32}, ArraySelectionParameter::AllowedComponentShapes{{1}})); @@ -67,10 +73,41 @@ Parameters ComputeAvgOrientationsFilter::parameters() const params.insert(std::make_unique(k_CellFeatureAttributeMatrixPath_Key, "Feature Attribute Matrix", "The path to the cell feature attribute matrix", DataPath({"Cell Feature Data"}))); params.insertSeparator(Parameters::Separator{"Output Feature Data"}); - params.insert(std::make_unique(k_AvgQuatsArrayName_Key, "Average Quaternions", - "The name of the array specifying the average orientation of the Feature in quaternion representation", "AvgQuats")); - params.insert(std::make_unique(k_AvgEulerAnglesArrayName_Key, "Average Euler Angles", - "The name of the array specifying the orientation of each Feature in Bunge convention (Z-X-Z)", "AvgEulerAngles")); + params.insert(std::make_unique(k_RodriguesQuatsArrayName_Key, "Average Rodrigues Quaternions", + "The name of the array specifying the average orientation based on the average Rodrigues vector of the Feature in quaternion representation", + "AvgQuats")); + params.insert(std::make_unique(k_RodriguesAvgEulerArrayName_Key, "Average Rodrigues Euler Angles", + "The name of the array specifying the orientation based on the average Rodrigues vector of each Feature in Bunge convention (Z-X-Z)", + "AvgEulerAngles")); + + params.insert(std::make_unique( + k_VonMisesFisherAvgQuatsArrayName_Key, "Average von Mises-Fisher Quaternions", + "The name of the array specifying the average orientation based on the von Mises-Fisher sampling of each Feature in quaternion representation", "vMF Avg Quats")); + params.insert(std::make_unique(k_VonMisesFisherAvgEulerArrayName_Key, "Average von Mises-Fisher Euler Angles", + "The name of the array specifying the average orientation based on the von Mises-Fisher sampling of each Feature in Bunge convention (Z-X-Z)", + "vMF Avg EulerAngles")); + params.insert(std::make_unique(k_VonMisesFisherKappaArrayName_Key, "von Mises-Fisher Kappa Values", + "The name of the array specifying the kappa values from the von Mises-Fisher sampling of each Feature", "vMF Kappas")); + + params.insert(std::make_unique(k_WatsonAvgQuatsArrayName_Key, "Average Watson Quaternions", + "The name of the array specifying the average orientation based on the Watson sampling of each Feature in quaternion representation", + "Watson Avg Quats")); + params.insert(std::make_unique(k_WatsonAvgEulerArrayName_Key, "Average Watson Angles", + "The name of the array specifying the average orientation based on the Watson sampling of each Feature in Bunge convention (Z-X-Z)", + "Watson Avg EulerAngles")); + params.insert(std::make_unique(k_WatsonKappaArrayName_Key, "Watson Kappa Values", + "The name of the array specifying the kappa values from the Watson sampling of each Feature", "Watson Kappas")); + + params.linkParameters(k_UseRodriguesAverage_Key, k_RodriguesAvgEulerArrayName_Key, true); + params.linkParameters(k_UseRodriguesAverage_Key, k_RodriguesQuatsArrayName_Key, true); + + params.linkParameters(k_UseVonMisesFisher_Key, k_VonMisesFisherAvgQuatsArrayName_Key, true); + params.linkParameters(k_UseVonMisesFisher_Key, k_VonMisesFisherAvgEulerArrayName_Key, true); + params.linkParameters(k_UseVonMisesFisher_Key, k_VonMisesFisherKappaArrayName_Key, true); + + params.linkParameters(k_UseWatson_Key, k_WatsonAvgQuatsArrayName_Key, true); + params.linkParameters(k_UseWatson_Key, k_WatsonAvgEulerArrayName_Key, true); + params.linkParameters(k_UseWatson_Key, k_WatsonKappaArrayName_Key, true); return params; } @@ -78,7 +115,9 @@ Parameters ComputeAvgOrientationsFilter::parameters() const //------------------------------------------------------------------------------ IFilter::VersionType ComputeAvgOrientationsFilter::parametersVersion() const { - return 1; + return 2; + // Version 2 adds the ability to compute the von Mises-Fisher average and the Watson sampling average + // Version 2 also adds the option to NOT compute the Eulers/Quats from the original algorithm } //------------------------------------------------------------------------------ @@ -96,8 +135,22 @@ IFilter::PreflightResult ComputeAvgOrientationsFilter::preflightImpl(const DataS auto pCellQuatsArrayPathValue = filterArgs.value(k_CellQuatsArrayPath_Key); auto pCrystalStructuresArrayPathValue = filterArgs.value(k_CrystalStructuresArrayPath_Key); auto pCellFeatureAttributeMatrixPathValue = filterArgs.value(k_CellFeatureAttributeMatrixPath_Key); - auto pAvgQuatsArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_AvgQuatsArrayName_Key)); - auto pAvgEulerAnglesArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_AvgEulerAnglesArrayName_Key)); + + auto pUseRodriguesAverage_Key = filterArgs.value(k_UseRodriguesAverage_Key); + auto pAvgQuatsArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_RodriguesQuatsArrayName_Key)); + auto pAvgEulerAnglesArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_RodriguesAvgEulerArrayName_Key)); + + auto pUseVonMisesFisher = filterArgs.value(k_UseVonMisesFisher_Key); + auto pVMFQuatsArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherAvgQuatsArrayName_Key)); + auto pVMFEulerAnglesArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherAvgEulerArrayName_Key)); + auto pVMFKappaArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherKappaArrayName_Key)); + + auto pUseWatson = filterArgs.value(k_UseWatson_Key); + auto pWatsonQuatsArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonAvgQuatsArrayName_Key)); + auto pWatsonEulerAnglesArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonAvgEulerArrayName_Key)); + auto pWatsonKappaArrayPathValue = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonKappaArrayName_Key)); + + Result resultOutputActions; std::vector dataPaths; @@ -115,16 +168,53 @@ IFilter::PreflightResult ComputeAvgOrientationsFilter::preflightImpl(const DataS const auto& cellFeatAttMatrix = dataStructure.getDataRefAs(pCellFeatureAttributeMatrixPathValue); // Create output DataStructure Items - auto tDims = cellFeatAttMatrix.getShape(); - auto createAvgQuatAction = std::make_unique(DataType::float32, tDims, std::vector{4}, pAvgQuatsArrayPathValue); - auto createAvgEulerAction = std::make_unique(DataType::float32, tDims, std::vector{3}, pAvgEulerAnglesArrayPathValue); + auto tupleShape = cellFeatAttMatrix.getShape(); + if(pUseRodriguesAverage_Key) + { + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{4}, pAvgQuatsArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{3}, pAvgEulerAnglesArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + } + + if(pUseVonMisesFisher) + { + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{4}, pVMFQuatsArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{3}, pVMFEulerAnglesArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{1}, pVMFKappaArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + } - OutputActions actions; - actions.appendAction(std::move(createAvgQuatAction)); - actions.appendAction(std::move(createAvgEulerAction)); + if(pUseWatson) + { + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{4}, pWatsonQuatsArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{3}, pWatsonEulerAnglesArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + { + auto createArrayAction = std::make_unique(DataType::float32, tupleShape, std::vector{1}, pWatsonKappaArrayPathValue); + resultOutputActions.value().appendAction(std::move(createArrayAction)); + } + } // Return both the resultOutputActions and the preflightUpdatedValues via std::move() - return {std::move(actions)}; + return {resultOutputActions}; } //------------------------------------------------------------------------------ @@ -138,8 +228,22 @@ Result<> ComputeAvgOrientationsFilter::executeImpl(DataStructure& dataStructure, inputValues.cellQuatsArrayPath = filterArgs.value(k_CellQuatsArrayPath_Key); inputValues.crystalStructuresArrayPath = filterArgs.value(k_CrystalStructuresArrayPath_Key); auto pCellFeatureAttributeMatrixPathValue = filterArgs.value(k_CellFeatureAttributeMatrixPath_Key); - inputValues.avgQuatsArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_AvgQuatsArrayName_Key)); - inputValues.avgEulerAnglesArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_AvgEulerAnglesArrayName_Key)); + + inputValues.useRodriguesAverage = filterArgs.value(k_UseRodriguesAverage_Key); + inputValues.avgQuatsArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_RodriguesQuatsArrayName_Key)); + inputValues.avgEulerAnglesArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_RodriguesAvgEulerArrayName_Key)); + + inputValues.useVonMisesAverage = filterArgs.value(k_UseVonMisesFisher_Key); + inputValues.VMFQuatsArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherAvgQuatsArrayName_Key)); + inputValues.VMFEulerAnglesArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherAvgEulerArrayName_Key)); + inputValues.VMFKappaArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_VonMisesFisherKappaArrayName_Key)); + + inputValues.useWatsonAverage = filterArgs.value(k_UseWatson_Key); + inputValues.WatsonQuatsArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonAvgQuatsArrayName_Key)); + inputValues.WatsonEulerAnglesArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonAvgEulerArrayName_Key)); + inputValues.WatsonKappaArrayPath = pCellFeatureAttributeMatrixPathValue.createChildPath(filterArgs.value(k_WatsonKappaArrayName_Key)); + + inputValues.RandomSeed = 43514; // Let the Algorithm instance do the work return ComputeAvgOrientations(dataStructure, messageHandler, shouldCancel, &inputValues)(); @@ -171,9 +275,10 @@ Result ComputeAvgOrientationsFilter::FromSIMPLJson(const nlohmann::js results.push_back(SIMPLConversion::ConvertParameter(args, json, SIMPL::k_QuatsArrayPathKey, k_CellQuatsArrayPath_Key)); results.push_back( SIMPLConversion::ConvertParameter(args, json, SIMPL::k_CrystalStructuresArrayPathKey, k_CrystalStructuresArrayPath_Key)); - results.push_back(SIMPLConversion::ConvertParameter(args, json, SIMPL::k_AvgQuatsArrayPathKey, k_AvgQuatsArrayName_Key)); results.push_back( - SIMPLConversion::ConvertParameter(args, json, SIMPL::k_AvgEulerAnglesArrayPathKey, k_AvgEulerAnglesArrayName_Key)); + SIMPLConversion::ConvertParameter(args, json, SIMPL::k_AvgQuatsArrayPathKey, k_RodriguesQuatsArrayName_Key)); + results.push_back( + SIMPLConversion::ConvertParameter(args, json, SIMPL::k_AvgEulerAnglesArrayPathKey, k_RodriguesAvgEulerArrayName_Key)); Result<> conversionResult = MergeResults(std::move(results)); diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.hpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.hpp index 509d2eec12..a66893f9c2 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.hpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeAvgOrientationsFilter.hpp @@ -28,10 +28,22 @@ class ORIENTATIONANALYSIS_EXPORT ComputeAvgOrientationsFilter : public IFilter static constexpr StringLiteral k_CellPhasesArrayPath_Key = "cell_phases_array_path"; static constexpr StringLiteral k_CellQuatsArrayPath_Key = "cell_quats_array_path"; static constexpr StringLiteral k_CrystalStructuresArrayPath_Key = "crystal_structures_array_path"; - static constexpr StringLiteral k_AvgQuatsArrayName_Key = "avg_quats_array_name"; - static constexpr StringLiteral k_AvgEulerAnglesArrayName_Key = "avg_euler_angles_array_name"; static constexpr StringLiteral k_CellFeatureAttributeMatrixPath_Key = "cell_feature_attribute_matrix_path"; + static constexpr StringLiteral k_UseRodriguesAverage_Key = "use_rodrigues_average"; + static constexpr StringLiteral k_RodriguesAvgEulerArrayName_Key = "avg_euler_angles_array_name"; + static constexpr StringLiteral k_RodriguesQuatsArrayName_Key = "avg_quats_array_name"; + + static constexpr StringLiteral k_UseVonMisesFisher_Key = "use_vmf_average"; + static constexpr StringLiteral k_VonMisesFisherAvgQuatsArrayName_Key = "vmf_quat_array_name"; + static constexpr StringLiteral k_VonMisesFisherAvgEulerArrayName_Key = "vmf_euler_array_name"; + static constexpr StringLiteral k_VonMisesFisherKappaArrayName_Key = "vmf_kappa_array_name"; + + static constexpr StringLiteral k_UseWatson_Key = "use_watson_average"; + static constexpr StringLiteral k_WatsonAvgQuatsArrayName_Key = "watson_quat_array_name"; + static constexpr StringLiteral k_WatsonAvgEulerArrayName_Key = "watson_euler_array_name"; + static constexpr StringLiteral k_WatsonKappaArrayName_Key = "watson_kappa_array_name"; + /** * @brief Reads SIMPL json and converts it simplnx Arguments. * @param json diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeFZQuaternionsFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeFZQuaternionsFilter.cpp index 20d3be405c..a77f6e6b2f 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeFZQuaternionsFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ComputeFZQuaternionsFilter.cpp @@ -13,165 +13,10 @@ #include "simplnx/Utilities/ParallelDataAlgorithm.hpp" #include "simplnx/Utilities/SIMPLConversion.hpp" -#include -#include #include using namespace nx::core; -namespace -{ -/** - * @brief The GenerateFZQuatsImpl class implements a threaded algorithm that computes the Fundamental Zone Quaternion - * for a given Quaternion and Laue Class (which is based from the crystalStructures array - */ -class GenerateMaskedFZQuatsImpl -{ -public: - GenerateMaskedFZQuatsImpl(const Float32AbstractDataStore& quats, const Int32AbstractDataStore& phases, const UInt32AbstractDataStore& crystalStructures, const int32 numPhases, - std::unique_ptr& goodVoxels, Float32AbstractDataStore& fzQuats, const std::atomic_bool& shouldCancel, std::atomic_int32_t& warningCount) - : m_Quats(quats) - , m_CellPhases(phases) - , m_CrystalStructures(crystalStructures) - , m_NumPhases(numPhases) - , m_GoodVoxels(goodVoxels) - , m_FZQuats(fzQuats) - , m_ShouldCancel(shouldCancel) - , m_WarningCount(warningCount) - { - } - - ~GenerateMaskedFZQuatsImpl() = default; - - /** - * @brief convert - * @param start - * @param end - */ - void convert(size_t start, size_t end) const - { - std::vector ops = ebsdlib::LaueOps::GetAllOrientationOps(); - int32 phase = 0; - size_t index = 0; - - for(size_t i = start; i < end; i++) - { - if(m_ShouldCancel) - { - break; - } - phase = m_CellPhases.getValue(i); - - // Sanity check the phase data to make sure we do not walk off the end of the array - if(phase >= m_NumPhases) - { - m_WarningCount++; - } - - // Output initialized to zero by default - index = i * 4; - if(phase < m_NumPhases && m_GoodVoxels->isTrue(i) && m_CrystalStructures.getValue(phase) < ebsdlib::CrystalStructure::LaueGroupEnd) - { - ebsdlib::QuatD quatD = ebsdlib::QuatD(m_Quats.getValue(index), m_Quats.getValue(index + 1), m_Quats.getValue(index + 2), m_Quats.getValue(index + 3)); // Makes a copy into q - auto xtal = static_cast(m_CrystalStructures.getValue(phase)); // get the Laue Group - quatD = ops[xtal]->getFZQuat(quatD); - m_FZQuats.setValue(index, static_cast(quatD.x())); - m_FZQuats.setValue(index + 1, static_cast(quatD.y())); - m_FZQuats.setValue(index + 2, static_cast(quatD.z())); - m_FZQuats.setValue(index + 3, static_cast(quatD.w())); - } - } - } - - void operator()(const Range& range) const - { - convert(range.min(), range.max()); - } - -private: - const Float32AbstractDataStore& m_Quats; - const Int32AbstractDataStore& m_CellPhases; - const UInt32AbstractDataStore& m_CrystalStructures; - const int32 m_NumPhases = 0; - std::unique_ptr& m_GoodVoxels; - Float32AbstractDataStore& m_FZQuats; - const std::atomic_bool& m_ShouldCancel; - std::atomic_int32_t& m_WarningCount; -}; - -class GenerateFZQuatsImpl -{ -public: - GenerateFZQuatsImpl(const Float32AbstractDataStore& quats, const Int32AbstractDataStore& phases, const UInt32AbstractDataStore& crystalStructures, const int32 numPhases, - Float32AbstractDataStore& fzQuats, const std::atomic_bool& shouldCancel, std::atomic_int32_t& warningCount) - : m_Quats(quats) - , m_CellPhases(phases) - , m_CrystalStructures(crystalStructures) - , m_NumPhases(numPhases) - , m_FZQuats(fzQuats) - , m_ShouldCancel(shouldCancel) - , m_WarningCount(warningCount) - { - } - - ~GenerateFZQuatsImpl() = default; - - /** - * @brief convert - * @param start - * @param end - */ - void convert(size_t start, size_t end) const - { - std::vector ops = ebsdlib::LaueOps::GetAllOrientationOps(); - int32 phase = 0; - size_t index = 0; - - for(size_t i = start; i < end; i++) - { - if(m_ShouldCancel) - { - break; - } - phase = m_CellPhases.getValue(i); - - // Sanity check the phase data to make sure we do not walk off the end of the array - if(phase >= m_NumPhases) - { - m_WarningCount++; - } - - // Output initialized to zero by default - index = i * 4; - if(phase < m_NumPhases && m_CrystalStructures.getValue(phase) < ebsdlib::CrystalStructure::LaueGroupEnd) - { - ebsdlib::QuatD quatD = ebsdlib::QuatD(m_Quats.getValue(index), m_Quats.getValue(index + 1), m_Quats.getValue(index + 2), m_Quats.getValue(index + 3)); // Makes a copy into q - auto xtal = static_cast(m_CrystalStructures.getValue(phase)); // get the Laue Group - quatD = ops[xtal]->getFZQuat(quatD); - m_FZQuats.setValue(index, static_cast(quatD.x())); - m_FZQuats.setValue(index + 1, static_cast(quatD.y())); - m_FZQuats.setValue(index + 2, static_cast(quatD.z())); - m_FZQuats.setValue(index + 3, static_cast(quatD.w())); - } - } - } - - void operator()(const Range& range) const - { - convert(range.min(), range.max()); - } - -private: - const Float32AbstractDataStore& m_Quats; - const Int32AbstractDataStore& m_CellPhases; - const UInt32AbstractDataStore& m_CrystalStructures; - const int32 m_NumPhases = 0; - Float32AbstractDataStore& m_FZQuats; - const std::atomic_bool& m_ShouldCancel; - std::atomic_int32_t& m_WarningCount; -}; -} // namespace - namespace nx::core { //------------------------------------------------------------------------------ @@ -289,76 +134,16 @@ IFilter::PreflightResult ComputeFZQuaternionsFilter::preflightImpl(const DataStr Result<> ComputeFZQuaternionsFilter::executeImpl(DataStructure& dataStructure, const Arguments& filterArgs, const PipelineFilter* pipelineNode, const MessageHandler& messageHandler, const std::atomic_bool& shouldCancel, const ExecutionContext& executionContext) const { - auto pUseGoodVoxelsValue = filterArgs.value(k_UseMask_Key); - auto pQuatsArrayPathValue = filterArgs.value(k_QuatsArrayPath_Key); - auto pCellPhasesArrayPathValue = filterArgs.value(k_CellPhasesArrayPath_Key); - auto pGoodVoxelsArrayPathValue = filterArgs.value(k_MaskArrayPath_Key); - auto pCrystalStructuresArrayPathValue = filterArgs.value(k_CrystalStructuresArrayPath_Key); - auto pFZQuatsArrayPathValue = pQuatsArrayPathValue.replaceName(filterArgs.value(k_FZQuatsArrayName_Key)); - - auto& phaseArray = dataStructure.getDataRefAs(pCellPhasesArrayPathValue); - auto& quatArray = dataStructure.getDataRefAs(pQuatsArrayPathValue); - auto& xtalArray = dataStructure.getDataRefAs(pCrystalStructuresArrayPathValue); - auto* maskArray = dataStructure.getDataAs(pGoodVoxelsArrayPathValue); - auto& fzQuatArray = dataStructure.getDataRefAs(pFZQuatsArrayPathValue); - std::atomic_int32_t warningCount = 0; - auto numPhases = static_cast(xtalArray.getNumberOfTuples()); - - typename IParallelAlgorithm::AlgorithmArrays algArrays; - algArrays.push_back(&phaseArray); - algArrays.push_back(&quatArray); - algArrays.push_back(&xtalArray); - algArrays.push_back(&fzQuatArray); - - if(pUseGoodVoxelsValue) - { - algArrays.push_back(maskArray); - } - - try - { - // Parallel algorithm - ParallelDataAlgorithm dataAlg; - dataAlg.setRange(0ULL, static_cast(quatArray.getNumberOfTuples())); - dataAlg.requireArraysInMemory(algArrays); - - if(pUseGoodVoxelsValue) - { - std::unique_ptr maskArrayPtr = nullptr; - try - { - maskArrayPtr = MaskCompareUtilities::InstantiateMaskCompare(dataStructure, pGoodVoxelsArrayPathValue); - } catch(const std::out_of_range& exception) - { - // This really should NOT be happening as the path was verified during preflight BUT we may be calling this from - // some other context that is NOT going through the normal nx::core::IFilter API of Preflight and Execute - return MakeErrorResult(-49003, fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", pGoodVoxelsArrayPathValue.toString())); - } - - dataAlg.execute(::GenerateMaskedFZQuatsImpl(quatArray.getDataStoreRef(), phaseArray.getDataStoreRef(), xtalArray.getDataStoreRef(), numPhases, maskArrayPtr, fzQuatArray.getDataStoreRef(), - shouldCancel, warningCount)); - } - else - { - dataAlg.execute( - ::GenerateFZQuatsImpl(quatArray.getDataStoreRef(), phaseArray.getDataStoreRef(), xtalArray.getDataStoreRef(), numPhases, fzQuatArray.getDataStoreRef(), shouldCancel, warningCount)); - } - - if(warningCount > 0) - { - std::string errorMessage = fmt::format("The Ensemble Phase information only references {} phase(s) but {} cell(s) had a phase value greater than {}. This indicates a problem with the input " - "cell phase data. DREAM3D-NX may have given INCORRECT RESULTS.", - numPhases - 1, warningCount.load(), numPhases - 1); - - return {MakeErrorResult<>(-49004, errorMessage)}; - } - } catch(const ebsdlib::method_not_implemented& e) - { - return {MakeErrorResult<>(-49005, fmt::format("EbsdLib threw an exception when computing the fundamental zone data. {}", e.what()))}; - } + ComputeFZQuaternionsInputValues inputValues; + inputValues.UseMask = filterArgs.value(k_UseMask_Key); + inputValues.InputQuatsArrayPath = filterArgs.value(k_QuatsArrayPath_Key); + inputValues.CellPhasesArrayPath = filterArgs.value(k_CellPhasesArrayPath_Key); + inputValues.CrystalStructuresArrayPath = filterArgs.value(k_CrystalStructuresArrayPath_Key); + inputValues.MaskArrayPath = filterArgs.value(k_MaskArrayPath_Key); + inputValues.OutputFzQuatsArrayName = filterArgs.value(k_FZQuatsArrayName_Key); - return {}; + return ComputeFZQuaternions(dataStructure, messageHandler, shouldCancel, &inputValues)(); } namespace diff --git a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ConvertOrientationsToVertexGeometryFilter.cpp b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ConvertOrientationsToVertexGeometryFilter.cpp index 58a40d3f5a..b80716df6e 100644 --- a/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ConvertOrientationsToVertexGeometryFilter.cpp +++ b/src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/ConvertOrientationsToVertexGeometryFilter.cpp @@ -7,8 +7,8 @@ #include "simplnx/DataStructure/Geometry/VertexGeom.hpp" #include "simplnx/Filter/Actions/CopyDataObjectAction.hpp" +#include "simplnx/Filter/Actions/CreateArrayAction.hpp" #include "simplnx/Filter/Actions/CreateVertexGeometryAction.hpp" -#include "simplnx/Filter/Actions/MoveDataAction.hpp" #include "simplnx/Parameters/ArraySelectionParameter.hpp" #include "simplnx/Parameters/BoolParameter.hpp" #include "simplnx/Parameters/ChoicesParameter.hpp" @@ -109,7 +109,7 @@ IFilter::PreflightResult ConvertOrientationsToVertexGeometryFilter::preflightImp { auto inputRepType = static_cast(filterArgs.value(k_InputType_Key)); auto inputOrientationsArrayPath = filterArgs.value(k_InputOrientationArrayPath_Key); - auto vertexPathsToCopy = filterArgs.value(k_CopyVertexPaths_Key); + auto sourceDataPaths = filterArgs.value(k_CopyVertexPaths_Key); auto convertToFundamentalZone = filterArgs.value(k_ConvertToFundamentalZone_Key); auto cellPhasesArrayPath = filterArgs.value(k_CellPhasesArrayPath_Key); auto crystalStructuresArrayPath = filterArgs.value(k_CrystalStructuresArrayPath_Key); @@ -145,33 +145,33 @@ IFilter::PreflightResult ConvertOrientationsToVertexGeometryFilter::preflightImp } // Create the Vertex Geometry - auto createVertexGeometryAction = - std::make_unique(outputVertexGeometryPath, inputOrientationsArray.getNumberOfTuples(), outputVertexAttrMatrixName, outputSharedVertexListName); - resultOutputActions.value().appendAction(std::move(createVertexGeometryAction)); + { + auto createVertexGeometryAction = + std::make_unique(outputVertexGeometryPath, inputOrientationsArray.getNumberOfTuples(), outputVertexAttrMatrixName, outputSharedVertexListName); + resultOutputActions.value().appendAction(std::move(createVertexGeometryAction)); + } DataPath vertexAttrMatrixPath = outputVertexGeometryPath.createChildPath(outputVertexAttrMatrixName); - for(const auto& vertexPathToCopy : vertexPathsToCopy) + for(const auto& sourceDataPath : sourceDataPaths) { - auto& vertexDataArray = dataStructure.getDataRefAs(vertexPathToCopy); - DataType type = vertexDataArray.getDataType(); - DataPath copyPath = vertexAttrMatrixPath.createChildPath(vertexDataArray.getName()); - auto numTuples = vertexDataArray.getNumberOfTuples(); - auto components = vertexDataArray.getNumberOfComponents(); - const std::string dataStoreFormat = vertexDataArray.getDataFormat(); + auto& sourceDataArray = dataStructure.getDataRefAs(sourceDataPath); + DataType type = sourceDataArray.getDataType(); + DataPath destinationDataPath = vertexAttrMatrixPath.createChildPath(sourceDataArray.getName()); + auto numTuples = sourceDataArray.getNumberOfTuples(); + auto components = sourceDataArray.getComponentShape(); + const std::string dataStoreFormat = sourceDataArray.getDataFormat(); if(numTuples != inputOrientationsArray.getNumberOfTuples()) { return {MakeErrorResult(-1004, fmt::format("Array at path {} only has {} tuples, but it MUST have {} tuples to be copied into output vertex attribute matrix at path {}!", - vertexPathToCopy.toString(), numTuples, inputOrientationsArray.getNumberOfTuples(), vertexAttrMatrixPath.toString()))}; + sourceDataPath.toString(), numTuples, inputOrientationsArray.getNumberOfTuples(), vertexAttrMatrixPath.toString()))}; } - auto action = std::make_unique(vertexPathToCopy, copyPath, std::vector{}); + // We are using the CreateArrayAction here instead of the CopyArrayAction because the + // CopyArrayAction will have the destination DataArray have the same tuple shape + // as the source array, but in this case we don't want that. So we use the CreateArrayAction + // instead and manually update the tuple shape later one. + auto action = std::make_unique(type, ShapeType{numTuples}, components, destinationDataPath, dataStoreFormat); resultOutputActions.value().appendAction(std::move(action)); - - // auto moveDataAction = std::make_unique(vertexPathToCopy, vertexAttrMatrixPath); - // resultOutputActions.value().appendAction(std::move(moveDataAction)); - - // auto action = std::make_unique(type, std::vector{numTuples}, std::vector{components}, copyPath, dataStoreFormat); - // resultOutputActions.value().appendAction(std::move(action)); } // Return both the resultOutputActions and the preflightUpdatedValues via std::move() @@ -186,7 +186,7 @@ Result<> ConvertOrientationsToVertexGeometryFilter::executeImpl(DataStructure& d inputValues.InputOrientationType = static_cast(filterArgs.value(k_InputType_Key)); inputValues.InputOrientationArrayPath = filterArgs.value(k_InputOrientationArrayPath_Key); - inputValues.CopyVertexArrayPaths = filterArgs.value(k_CopyVertexPaths_Key); + inputValues.DataPathCopySources = filterArgs.value(k_CopyVertexPaths_Key); inputValues.ConvertToFundamentalZone = filterArgs.value(k_ConvertToFundamentalZone_Key); inputValues.CellPhasesArrayPath = filterArgs.value(k_CellPhasesArrayPath_Key); inputValues.CrystalStructuresArrayPath = filterArgs.value(k_CrystalStructuresArrayPath_Key); diff --git a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt index 5ba060ddbd..41a6a49f7c 100644 --- a/src/Plugins/OrientationAnalysis/test/CMakeLists.txt +++ b/src/Plugins/OrientationAnalysis/test/CMakeLists.txt @@ -130,7 +130,7 @@ if(EXISTS "${DREAM3D_DATA_DIR}" AND SIMPLNX_DOWNLOAD_TEST_FILES) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 6_6_stats_test_v2.tar.gz SHA512 e84999dec914d81efce4fc4237c49c9bf32e48381b1e79f58aa4df934f0d7606cd7a948f9a5e7b17a126a7944cc531b531cfdc70756ca3e2207b20734e089723) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_2_AvgCAxis.tar.gz SHA512 054d1fbb92baacfa79f0bd326ae32b5bfc67d0935413ea8f194980527ff9694bd21a49f73e43d59b731a567ea89190523176bf5c98e8585e7b206265d5b05143) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_compute_triangle_shapes_test.tar.gz SHA512 562ec65d22771817655ac85e267b1ea9822ab478ef732a52ce05052e4a171fc131bb879c9982df32b726d68a350f092695af97633f69e97226a6147cd0608b82) - download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_ComputeAvgOrientation.tar.gz SHA512 5ee425b9a45da4daab1d43a870b54d57ac8d8e300acda162fa7a974a22b9f5662a98e768cebac1844eb87fcb2bb01624ffaf117d3aacc722aa8701712cbf3c52) + download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_ComputeAvgOrientation_v2.tar.gz SHA512 2c2a691f1da301c449c20bafec65512d5134db38384ac7cb4c910880ccd87a260a5f011e905f35b97abff3952309f109c737c63ec3c833708926827a62a92efc) download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME 7_read_oem_ebsd_h5_files.tar.gz SHA512 95fa4cc0fb6bce26bfd6d28c0205a9e0f6497ad876a7cc4381117e668836936b35683ed5ec757cadd124a6c6fec7b38fe11f72512105bbd677bcf4210892ac19) 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 align_sections_misorientation.tar.gz SHA512 8a186b2e96dd94a8583eacaec768c252885d89c8f5734b6511d573235beae075971e6e81b42bb517b7cd617fc478ed394abf8ea4fe3188f50d340f90573013f4) diff --git a/src/Plugins/OrientationAnalysis/test/ComputeAvgOrientationsTest.cpp b/src/Plugins/OrientationAnalysis/test/ComputeAvgOrientationsTest.cpp index 43838cb841..ab9251379b 100644 --- a/src/Plugins/OrientationAnalysis/test/ComputeAvgOrientationsTest.cpp +++ b/src/Plugins/OrientationAnalysis/test/ComputeAvgOrientationsTest.cpp @@ -60,19 +60,32 @@ const DataPath k_FeatureAttrMatPath = k_ImageDataContainerPath.createChildPath(" const std::string k_AvgQuatsName = "Computed AvgQuats"; const std::string k_AvgEulersName = "Computed AvgEulerAngles"; +const std::string k_ComputedWatsonAvgQuatsName = "Computed Watson AvgQuats"; +const std::string k_ComputedWatsonAvgEulersName = "Computed Watson AvgEulerAngles"; +const std::string k_ComputedWatsonKappasName = "Computed Watson Kappas"; + +const std::string k_ComputedVMFAvgQuatsName = "Computed VMF AvgQuats"; +const std::string k_ComputedVMFAvgEulersName = "Computed VMF AvgEulerAngles"; +const std::string k_ComputedVMFKappasName = "Computed VMF Kappas"; + const std::string k_Exemplar_AvgQuatsName = "AvgQuats"; const std::string k_Exemplar_AvgEulersName = "AvgEulerAngles"; +const std::string k_Exemplar_WatsonAvgQuatsName = "Watson Avg Quats"; +const std::string k_Exemplar_WatsonAvgEulersName = "Watson Avg EulerAngles"; + +const std::string k_Exemplar_VMFAvgQuatsName = "vMF Avg Quats"; +const std::string k_Exemplar_VMFAvgEulersName = "vMF Avg EulerAngles"; } // namespace compute_avg_orientation TEST_CASE("OrientationAnalysis::ComputeAvgOrientations", "[OrientationAnalysis][ComputeAvgOrientationsFilter]") { UnitTest::LoadPlugins(); - const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "7_ComputeAvgOrientation.tar.gz", "7_ComputeAvgOrientation"); + const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "7_ComputeAvgOrientation_v2.tar.gz", "7_ComputeAvgOrientation_v2"); // Read Exemplar DREAM3D File Filter - auto exemplarFilePath = fs::path(fmt::format("{}/7_ComputeAvgOrientation/7_ComputeAvgOrientation.dream3d", unit_test::k_TestFilesDir)); + auto exemplarFilePath = fs::path(fmt::format("{}/7_ComputeAvgOrientation_v2/7_ComputeAvgOrientation_v2.dream3d", unit_test::k_TestFilesDir)); DataStructure dataStructure = UnitTest::LoadDataStructure(exemplarFilePath); // Instantiate the filter, a DataStructure object and an Arguments Object @@ -84,10 +97,22 @@ TEST_CASE("OrientationAnalysis::ComputeAvgOrientations", "[OrientationAnalysis][ args.insertOrAssign(ComputeAvgOrientationsFilter::k_CellPhasesArrayPath_Key, std::make_any(compute_avg_orientation::k_PhasesPath)); args.insertOrAssign(ComputeAvgOrientationsFilter::k_CellQuatsArrayPath_Key, std::make_any(compute_avg_orientation::k_QuatsPath)); args.insertOrAssign(ComputeAvgOrientationsFilter::k_CrystalStructuresArrayPath_Key, std::make_any(compute_avg_orientation::k_CrystalStructuresPath)); - args.insertOrAssign(ComputeAvgOrientationsFilter::k_AvgQuatsArrayName_Key, std::make_any(compute_avg_orientation::k_AvgQuatsName)); - args.insertOrAssign(ComputeAvgOrientationsFilter::k_AvgEulerAnglesArrayName_Key, std::make_any(compute_avg_orientation::k_AvgEulersName)); args.insertOrAssign(ComputeAvgOrientationsFilter::k_CellFeatureAttributeMatrixPath_Key, std::make_any(compute_avg_orientation::k_FeatureAttrMatPath)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_UseRodriguesAverage_Key, std::make_any(true)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_RodriguesQuatsArrayName_Key, std::make_any(compute_avg_orientation::k_AvgQuatsName)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_RodriguesAvgEulerArrayName_Key, std::make_any(compute_avg_orientation::k_AvgEulersName)); + + args.insertOrAssign(ComputeAvgOrientationsFilter::k_UseWatson_Key, std::make_any(true)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_WatsonAvgQuatsArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedWatsonAvgQuatsName)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_WatsonAvgEulerArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedWatsonAvgEulersName)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_WatsonKappaArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedWatsonKappasName)); + + args.insertOrAssign(ComputeAvgOrientationsFilter::k_UseVonMisesFisher_Key, std::make_any(true)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_VonMisesFisherAvgQuatsArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedVMFAvgQuatsName)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_VonMisesFisherAvgEulerArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedVMFAvgEulersName)); + args.insertOrAssign(ComputeAvgOrientationsFilter::k_VonMisesFisherKappaArrayName_Key, std::make_any(compute_avg_orientation::k_ComputedVMFKappasName)); + // Preflight the filter and check result auto preflightResult = filter.preflight(dataStructure, args); SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions); @@ -100,15 +125,41 @@ TEST_CASE("OrientationAnalysis::ComputeAvgOrientations", "[OrientationAnalysis][ WriteTestDataStructure(dataStructure, fmt::format("{}/compute_average_orientations.dream3d", unit_test::k_BinaryTestOutputDir)); #endif - DataPath computedAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_AvgEulersName); - DataPath exemplarAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_AvgEulersName); + { + DataPath computedAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_AvgEulersName); + DataPath exemplarAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_AvgEulersName); + + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgEulersPath, exemplarAvgEulersPath, 5.0E-7f, false); + + DataPath computedAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_AvgQuatsName); + DataPath exemplarAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_AvgQuatsName); + + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgQuatsPath, exemplarAvgQuatsPath, 5.0E-7f, false); + } + + { + DataPath computedAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_ComputedWatsonAvgEulersName); + DataPath exemplarAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_WatsonAvgEulersName); + + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgEulersPath, exemplarAvgEulersPath, 5.0E-7f, false); + + DataPath computedAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_ComputedWatsonAvgQuatsName); + DataPath exemplarAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_WatsonAvgQuatsName); + + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgQuatsPath, exemplarAvgQuatsPath, 5.0E-7f, false); + } + + { + DataPath computedAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_ComputedVMFAvgEulersName); + DataPath exemplarAvgEulersPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_VMFAvgEulersName); - UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgEulersPath, exemplarAvgEulersPath, 5.0E-7f, false); + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgEulersPath, exemplarAvgEulersPath, 5.0E-7f, false); - DataPath computedAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_AvgQuatsName); - DataPath exemplarAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_AvgQuatsName); + DataPath computedAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_ComputedVMFAvgQuatsName); + DataPath exemplarAvgQuatsPath = compute_avg_orientation::k_FeatureAttrMatPath.createChildPath(compute_avg_orientation::k_Exemplar_VMFAvgQuatsName); - UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgQuatsPath, exemplarAvgQuatsPath, 5.0E-7f, false); + UnitTest::CompareFloatArraysWithNans(dataStructure, computedAvgQuatsPath, exemplarAvgQuatsPath, 5.0E-7f, false); + } UnitTest::CheckArraysInheritTupleDims(dataStructure); } diff --git a/test/ConversionTest/BackwardsCompatibilityTest.cpp b/test/ConversionTest/BackwardsCompatibilityTest.cpp index 747820a133..e21468e9b6 100644 --- a/test/ConversionTest/BackwardsCompatibilityTest.cpp +++ b/test/ConversionTest/BackwardsCompatibilityTest.cpp @@ -725,7 +725,14 @@ const std::map> k_KeyIgnoreMap = { "max_suffix_name", "mean_suffix_name", "min_suffix_name", "std_deviation_suffix_name", "summation_suffix_name"}}, // ExtractVertexGeometryFilter std::pair>{Uuid::FromString("621a71ca-124b-4471-ad1a-02f05ffba099").value(), - std::vector{"output_shared_vertex_list_name", "output_vertex_attr_matrix_name"}}}; + std::vector{"output_shared_vertex_list_name", "output_vertex_attr_matrix_name"}}, + + // ComputeAvgOrientationsFilter + std::pair>{Uuid::FromString("086ddb9a-928f-46ab-bad6-b1498270d71e").value(), + std::vector{"use_vmf_average", "use_watson_average", "vmf_euler_array_name", "vmf_kappa_array_name", "vmf_quat_array_name", + "watson_euler_array_name", "watson_kappa_array_name", "watson_quat_array_name"}} + +}; } // namespace diff --git a/test/UnitTestCommon/include/simplnx/UnitTest/UnitTestCommon.hpp b/test/UnitTestCommon/include/simplnx/UnitTest/UnitTestCommon.hpp index f7954f4c4b..67572b5159 100644 --- a/test/UnitTestCommon/include/simplnx/UnitTest/UnitTestCommon.hpp +++ b/test/UnitTestCommon/include/simplnx/UnitTest/UnitTestCommon.hpp @@ -1663,6 +1663,7 @@ inline void CheckArraysInheritTupleDims(const DataStructure& dataStructure, cons REQUIRE(daPathsOpt.has_value()); for(const auto& daPath : daPathsOpt.value()) { + INFO(fmt::format("AttributeMatrix: {} Array: '{}'", amPath.toString(), daPath.toString())); const auto& arr = dataStructure.getDataAs(daPath); REQUIRE(attrMatrix.getShape() == arr->getTupleShape()); }