Skip to content

ENH: Updates codes to use Eigen instead of built in Matrix Operations #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: develop
Choose a base branch
from
Open

ENH: Updates codes to use Eigen instead of built in Matrix Operations #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 42 additions & 49 deletions src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,14 @@

#include "simplnx/Common/Constants.hpp"
#include "simplnx/DataStructure/DataArray.hpp"
#include "simplnx/DataStructure/DataGroup.hpp"
#include "simplnx/DataStructure/NeighborList.hpp"
#include "simplnx/Utilities/Math/GeometryMath.hpp"
#include "simplnx/Utilities/Math/MatrixMath.hpp"

#include "EbsdLib/LaueOps/LaueOps.h"

#include "Plugins/OrientationAnalysis/src/OrientationAnalysis/utilities/OrientationUtilities.hpp"

#include <random>

using namespace nx::core;
Expand All @@ -19,15 +20,15 @@ void RandomizeFeatureIds(usize totalPoints, usize totalFeatures, Int32Array& cel
{
// Generate an even distribution of numbers between the min and max range
std::mt19937_64 gen(seed);
std::uniform_int_distribution<int64> dist(0, totalFeatures - 1);
std::uniform_int_distribution<int64> dist(0, static_cast<int64>(totalFeatures) - 1);

std::vector<int32> 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<int32>(dist(gen)); // Random remaining position.
const auto r = static_cast<int32>(dist(gen)); // Random remaining position.
if(r >= totalFeatures)
{
continue;
Expand Down Expand Up @@ -94,10 +95,10 @@ void GroupMicroTextureRegions::execute()
int32 seed = 0;
int32 list1size = 0, list2size = 0, listsize = 0;
int32 neigh = 0;
bool patchGrouping = false;

while(seed >= 0)
{
bool patchGrouping = false;
parentcount++;
seed = getSeed(parentcount);
if(seed >= 0)
Expand All @@ -106,7 +107,7 @@ void GroupMicroTextureRegions::execute()
for(std::vector<int32>::size_type j = 0; j < grouplist.size(); j++)
{
int32 firstfeature = grouplist[j];
list1size = int32(neighborlist[firstfeature].size());
list1size = static_cast<int32>(neighborlist[firstfeature].size());
if(m_InputValues->UseNonContiguousNeighbors)
{
list2size = nonContigNeighList->getListSize(firstfeature);
Expand Down Expand Up @@ -156,7 +157,7 @@ void GroupMicroTextureRegions::execute()
for(std::vector<int32_t>::size_type j = 0; j < grouplist.size(); j++)
{
int32_t firstfeature = grouplist[j];
listsize = int32_t(neighborlist[firstfeature].size());
listsize = static_cast<int32_t>(neighborlist[firstfeature].size());
for(int32_t l = 0; l < listsize; l++)
{
neigh = neighborlist[firstfeature][l];
Expand All @@ -179,12 +180,10 @@ void GroupMicroTextureRegions::execute()
// -----------------------------------------------------------------------------
Result<> GroupMicroTextureRegions::operator()()
{
// m_AvgCAxes = {0.0f, 0.0f, 0.0f};
m_Generator = std::mt19937_64(m_InputValues->SeedValue);
m_Distribution = std::uniform_real_distribution<float32>(0.0f, 1.0f);

m_AvgCAxes[0] = 0.0f;
m_AvgCAxes[1] = 0.0f;
m_AvgCAxes[2] = 0.0f;
auto& featureParentIds = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeatureParentIdsArrayName);
featureParentIds.fill(-1);

Expand All @@ -205,7 +204,7 @@ Result<> GroupMicroTextureRegions::operator()()
cellParentIds[k] = featureParentIds[featureIds[k]];
}

// By default we randomize grains !!! COMMENT OUT FOR DEMONSTRATION !!!
// By default, we randomize grains !!! COMMENT OUT FOR DEMONSTRATION !!!
m_MessageHandler(IFilter::Message::Type::Info, "Randomizing Parent Ids");
RandomizeFeatureIds(totalPoints, m_NumTuples, cellParentIds, featureParentIds, featureIds, m_InputValues->SeedValue);

Expand All @@ -219,9 +218,6 @@ int GroupMicroTextureRegions::getSeed(int32 newFid)
auto& featurePhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeaturePhasesArrayPath);

usize numFeatures = featurePhases.getNumberOfTuples();

float32 g1[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
float32 g1t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
int32 voxelSeed = -1;

// Precalculate some constants
Expand Down Expand Up @@ -250,6 +246,7 @@ int GroupMicroTextureRegions::getSeed(int32 newFid)
// std::ofstream fout ("/tmp/GroupMicroTextureInitialVoxelSeeds.txt", std::ios_base::out | std::ios_base::app);
// fout << fmt::format("Feature Parent Id: {} | X: {}, Y: {}\n", voxelSeed, centroids.getComponent(voxelSeed, 0), centroids.getComponent(voxelSeed, 1));
// }
const Eigen::Vector3f cAxis = {0.0f, 0.0f, 1.0f};

if(voxelSeed >= 0)
{
Expand All @@ -261,20 +258,22 @@ int GroupMicroTextureRegions::getSeed(int32 newFid)
auto& volumes = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->VolumesArrayPath);
auto& avgQuats = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AvgQuatsArrayPath);

usize index = voxelSeed * 4;
OrientationTransformation::qu2om<QuatF, OrientationF>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g1);

std::array<float32, 3> c1 = {0.0f, 0.0f, 0.0f};
std::array<float32, 3> cAxis = {0.0f, 0.0f, 1.0f};
// transpose the g matrix so when c-axis is multiplied by it,
const usize index = voxelSeed * 4;
const auto om = OrientationTransformation::qu2om<QuatF, OrientationF>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]});
// transpose the g matrix so when it multiplies c-axis,
// it will give the sample direction that the c-axis is along
MatrixMath::Transpose3x3(g1, g1t);
MatrixMath::Multiply3x3with3x1(g1t, cAxis.data(), c1.data());
const OrientationUtilities::Matrix3fR g1t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om);

Eigen::Vector3f c1 = g1t * cAxis;
// normalize so that the dot product can be taken below without
// dividing by the magnitudes (they would be 1)
MatrixMath::Normalize3x1(c1.data());
MatrixMath::Copy3x1(c1.data(), m_AvgCAxes.data());
MatrixMath::Multiply3x1withConstant(m_AvgCAxes.data(), volumes.getValue(voxelSeed));
c1.normalize();

c1 = c1 * volumes[voxelSeed];

m_RunningAvgOrientation[0] = c1[0];
m_RunningAvgOrientation[1] = c1[1];
m_RunningAvgOrientation[2] = c1[2];
}
}

Expand All @@ -284,57 +283,51 @@ int GroupMicroTextureRegions::getSeed(int32 newFid)
// -----------------------------------------------------------------------------
bool GroupMicroTextureRegions::determineGrouping(int32 referenceFeature, int32 neighborFeature, int32 newFid)
{
uint32 phase1 = 0;
float32 g1[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
float32 g2[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
float32 g1t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
float32 g2t[3][3] = {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}};
std::array<float32, 3> c1 = {0.0f, 0.0f, 0.0f};
std::array<float32, 3> caxis = {0.0f, 0.0f, 1.0f};
Eigen::Vector3f cAxis = {0.0f, 0.0f, 1.0f};

auto& featurePhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeaturePhasesArrayPath);
auto& featureParentIds = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeatureParentIdsArrayName);
auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
if(featureParentIds[neighborFeature] == -1 && featurePhases[referenceFeature] > 0 && featurePhases[neighborFeature] > 0)
{
uint32 phase1 = 0;
Eigen::Vector3f c1;
auto& avgQuats = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AvgQuatsArrayPath);
if(!m_InputValues->UseRunningAverage)
{
usize index = referenceFeature * 4;
phase1 = crystalStructures[featurePhases[referenceFeature]];
OrientationTransformation::qu2om<QuatF, Orientation<float32>>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g1);

// transpose the g matrix so when c-axis is multiplied by it,
const auto om = OrientationTransformation::qu2om<QuatF, OrientationF>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]});
// transpose the g matrix so when it multiplies c-axis,
// it will give the sample direction that the c-axis is along
MatrixMath::Transpose3x3(g1, g1t);
MatrixMath::Multiply3x3with3x1(g1t, caxis.data(), c1.data());
const OrientationUtilities::Matrix3fR g1t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om);
c1 = g1t * cAxis;
// normalize so that the dot product can be taken below without
// dividing by the magnitudes (they would be 1)
MatrixMath::Normalize3x1(c1.data());
c1.normalize();
}
uint32 phase2 = crystalStructures[featurePhases[neighborFeature]];
if(phase1 == phase2 && (phase1 == EbsdLib::CrystalStructure::Hexagonal_High))
{
usize index = neighborFeature * 4;
OrientationTransformation::qu2om<QuatF, OrientationF>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]}).toGMatrix(g2);

std::array<float32, 3> c2 = {0.0f, 0.0f, 0.0f};
// transpose the g matrix so when c-axis is multiplied by it,
const auto om = OrientationTransformation::qu2om<QuatF, OrientationF>({avgQuats[index + 0], avgQuats[index + 1], avgQuats[index + 2], avgQuats[index + 3]});
// transpose the g matrix so when it multiplies c-axis,
// it will give the sample direction that the c-axis is along
MatrixMath::Transpose3x3(g2, g2t);
MatrixMath::Multiply3x3with3x1(g2t, caxis.data(), c2.data());
const OrientationUtilities::Matrix3fR g2t = OrientationUtilities::OrientationMatrixToGMatrixTranspose(om);

Eigen::Vector3f c2 = g2t * cAxis;
// normalize so that the dot product can be taken below without
// dividing by the magnitudes (they would be 1)
MatrixMath::Normalize3x1(c2.data());
c2.normalize();

float32 w;
if(m_InputValues->UseRunningAverage)
{
w = GeometryMath::CosThetaBetweenVectors(Point3Df{m_AvgCAxes}, Point3Df{c2});
w = GeometryMath::CosThetaBetweenVectors(m_RunningAvgOrientation, c2);
}
else
{
w = GeometryMath::CosThetaBetweenVectors(Point3Df{c1}, Point3Df{c2});
w = GeometryMath::CosThetaBetweenVectors(c1, c2);
}

if(w < -1.0f)
Expand All @@ -348,15 +341,15 @@ bool GroupMicroTextureRegions::determineGrouping(int32 referenceFeature, int32 n
w = std::acos(w);

// Convert user defined tolerance to radians.
float32 cAxisToleranceRad = m_InputValues->CAxisTolerance * nx::core::Constants::k_PiD / 180.0f;
float32 cAxisToleranceRad = m_InputValues->CAxisTolerance * nx::core::Constants::k_PiF / 180.0f;
if(w <= cAxisToleranceRad || (nx::core::Constants::k_PiD - w) <= cAxisToleranceRad)
{
featureParentIds[neighborFeature] = newFid;
if(m_InputValues->UseRunningAverage)
{
auto& volumes = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->VolumesArrayPath);
MatrixMath::Multiply3x1withConstant(c2.data(), volumes.getValue(neighborFeature));
MatrixMath::Add3x1s(m_AvgCAxes.data(), c2.data(), m_AvgCAxes.data());
MatrixMath::Multiply3x1withConstant(c2.data(), volumes[neighborFeature]);
MatrixMath::Add3x1s(m_RunningAvgOrientation.data(), c2.data(), m_RunningAvgOrientation.data());
}
return true;
}
Expand Down
8 changes: 3 additions & 5 deletions src/SimplnxReview/Filters/Algorithms/GroupMicroTextureRegions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,8 @@
#include "simplnx/DataStructure/DataPath.hpp"
#include "simplnx/DataStructure/DataStructure.hpp"
#include "simplnx/Filter/IFilter.hpp"
#include "simplnx/Parameters/ArraySelectionParameter.hpp"
#include "simplnx/Parameters/BoolParameter.hpp"
#include "simplnx/Parameters/NumberParameter.hpp"
#include "simplnx/Parameters/StringParameter.hpp"

#include <Eigen/Dense>

#include <random>

Expand Down Expand Up @@ -67,7 +65,7 @@ class SIMPLNXREVIEW_EXPORT GroupMicroTextureRegions
const IFilter::MessageHandler& m_MessageHandler;

usize m_NumTuples = 0;
std::array<float32, 3> m_AvgCAxes = {0.0f, 0.0f, 0.0f};
Eigen::Vector3<float32> m_RunningAvgOrientation = {0.0f, 0.0f, 0.0f};
std::mt19937_64 m_Generator = {};
std::uniform_real_distribution<float32> m_Distribution = {};
};
Expand Down