Repository navigation
BUG: Compute C-Axis Locations writes a zero vector for a zero quaternion #1891
Description
Problem
Compute C-Axis Locations (ComputeCAxisLocationsFilter) writes the vector [0, 0, 0] for a tuple whose quaternion is all zeros and whose phase is hexagonal. A zero quaternion defines no rotation, so it has no c-axis direction. The expected output is NaN in all three components, the same value that the filter writes for non-hexagonal phases. The output [0, 0, 0] is not a unit vector, but later filters and plots can use it as a valid direction without an error.
The defect occurs at commit 705a1b7 on develop.
Steps to Reproduce
- Make a Float32 quaternion array with 2 tuples and 4 components: tuple 0 is
[0, 0, 0, 1](identity), tuple 1 is[0, 0, 0, 0]. - Make a UInt32 Crystal Structures array with 2 tuples:
[UnknownCrystalStructure, Hexagonal_High]. - Make an Int32 Phases array with 2 tuples, both set to
1. - Run Compute C-Axis Locations with these three arrays.
- Read the output array. Tuple 0 is
[0, 0, 1], which is correct. Tuple 1 is[0, 0, 0]. It should be[NaN, NaN, NaN].
Cause
For hexagonal tuples, the algorithm rotates the crystal [001] axis into the sample frame and normalizes it:
const usize qi = chunkTupleIdx * 4;
ebsdlib::OrientationMatrixFType oMatrix = ebsdlib::QuaternionFType(quatBuf[qi], quatBuf[qi + 1], quatBuf[qi + 2], quatBuf[qi + 3]).toOrientationMatrix();
c1 = oMatrix.transpose() * cAxis;
c1.normalize();
if(c1[2] < 0)
{
c1 *= -1.0f;
}
outputBuf[chunkTupleIdx * 3] = c1[0];For a zero quaternion, every element of the orientation matrix is zero, so c1 is the zero vector. Eigen's normalize() does not change a vector whose norm is zero, so [0, 0, 0] goes to the output. The algorithm does not check the norm before it writes the result.
Verification Unit Test
Add this test to src/Plugins/OrientationAnalysis/test/ComputeCAxisLocationsTest.cpp:
TEST_CASE("OrientationAnalysis::ComputeCAxisLocationsFilter: Zero quaternion gives NaN", "[OrientationAnalysis][ComputeCAxisLocationsFilter]")
{
UnitTest::LoadPlugins();
DataStructure dataStructure;
// Tuple 0 is the identity rotation. Tuple 1 is the all-zero quaternion, which defines no rotation.
Float32Array* quats = UnitTest::CreateTestDataArray<float32>(dataStructure, "Quats", {2}, {4});
const std::array<float32, 8> quatValues = {0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f};
for(usize i = 0; i < quatValues.size(); i++)
{
quats->setValue(i, quatValues[i]);
}
UInt32Array* crystalStructures = UnitTest::CreateTestDataArray<uint32>(dataStructure, "CrystalStructures", {2}, {1});
crystalStructures->setValue(0, ebsdlib::CrystalStructure::UnknownCrystalStructure);
crystalStructures->setValue(1, ebsdlib::CrystalStructure::Hexagonal_High);
Int32Array* phases = UnitTest::CreateTestDataArray<int32>(dataStructure, "Phases", {2}, {1});
phases->setValue(0, 1);
phases->setValue(1, 1);
DataPath cAxisLocationsPath({"CAxisLocations"});
ComputeCAxisLocationsFilter filter;
Arguments args;
args.insertOrAssign(ComputeCAxisLocationsFilter::k_QuatsArrayPath_Key, std::make_any<DataPath>(std::vector<std::string>{"Quats"}));
args.insertOrAssign(ComputeCAxisLocationsFilter::k_CellPhasesArrayPath_Key, std::make_any<DataPath>(std::vector<std::string>{"Phases"}));
args.insertOrAssign(ComputeCAxisLocationsFilter::k_CrystalStructuresArrayPath_Key, std::make_any<DataPath>(std::vector<std::string>{"CrystalStructures"}));
args.insertOrAssign(ComputeCAxisLocationsFilter::k_CAxisLocationsArrayName_Key, std::make_any<std::string>(cAxisLocationsPath.getTargetName()));
auto preflightResult = filter.preflight(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions);
auto executeResult = filter.execute(dataStructure, args);
SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result);
REQUIRE_NOTHROW(dataStructure.getDataRefAs<Float32Array>(cAxisLocationsPath));
auto& cAxisLocations = dataStructure.getDataRefAs<Float32Array>(cAxisLocationsPath);
// The identity rotation keeps the c-axis on +Z.
REQUIRE(cAxisLocations.getComponent(0, 0) == Approx(0.0f).margin(1e-7f));
REQUIRE(cAxisLocations.getComponent(0, 1) == Approx(0.0f).margin(1e-7f));
REQUIRE(cAxisLocations.getComponent(0, 2) == Approx(1.0f).margin(1e-7f));
// A zero quaternion has no c-axis direction, so all three components must be NaN,
// the same value the filter writes for non-hexagonal phases.
for(usize j = 0; j < 3; j++)
{
INFO(fmt::format("component {} = {}", j, cAxisLocations.getComponent(1, j)));
REQUIRE(std::isnan(cAxisLocations.getComponent(1, j)));
}
UnitTest::CheckArraysInheritTupleDims(dataStructure);
}Run only this test:
ctest -R "ComputeCAxisLocationsFilter: Zero quaternion gives NaN" --output-on-failure
Failure output (Ubuntu 22.04, gcc-13, ci-linux-x64 preset, commit 705a1b7):
src/Plugins/OrientationAnalysis/test/ComputeCAxisLocationsTest.cpp:543: FAILED:
REQUIRE( std::isnan(cAxisLocations.getComponent(1, j)) )
with expansion:
false
with message:
component 0 = 0
===============================================================================
test cases: 1 | 1 failed
assertions: 8 | 7 passed | 1 failed
The other 7 ComputeCAxisLocationsFilter tests pass in the same run.
Proposed Fix
In ComputeCAxisLocations::operator()() (src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeCAxisLocations.cpp), check the length of c1 before it is normalized. If the length is zero or not finite, write NaN for the three components, as the else branch does for non-hexagonal phases:
c1 = oMatrix.transpose() * cAxis;
- c1.normalize();
+ const float32 cAxisNorm = c1.norm();
+ if(!std::isfinite(cAxisNorm) || cAxisNorm == 0.0f)
+ {
+ outputBuf[chunkTupleIdx * 3] = NAN;
+ outputBuf[chunkTupleIdx * 3 + 1] = NAN;
+ outputBuf[chunkTupleIdx * 3 + 2] = NAN;
+ continue;
+ }
+ c1 /= cAxisNorm;
if(c1[2] < 0)Update the filter documentation to state that a zero or non-finite quaternion gives NaN.
Notes
- This defect was found during the review of FILT: Add Compute C-Axis Reference Misalignments filter #1890. The new Compute C-Axis Reference Misalignments filter in that pull request had the same problem, and the pull request fixes it there with the same norm check.
- BUG: C16: Degenerate input not handled in 39 filters #1875 lists degenerate-input defects in other filters, but its body does not list this filter.
- A quaternion with a NaN or infinite component probably also gives a wrong result in this filter. That case was not tested for this issue.
Activity