diff --git a/Source/EbsdLib/Core/DirectionalStats.cpp b/Source/EbsdLib/Core/DirectionalStats.cpp index ff6167ed..ccd1a365 100644 --- a/Source/EbsdLib/Core/DirectionalStats.cpp +++ b/Source/EbsdLib/Core/DirectionalStats.cpp @@ -1,3 +1,101 @@ +/* ============================================================================ + * Copyright (c) 2026 BlueQuartz Software, LLC + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * Redistributions in binary form must reproduce the above copyright notice, this + * list of conditions and the following disclaimer in the documentation and/or + * other materials provided with the distribution. + * + * Neither the name of BlueQuartz Software, the US Air Force, nor the names of its + * contributors may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + +/* ============================================================================ + * THIRD-PARTY ATTRIBUTION: EMsoft + * + * This file is a C++ port of the directional statistics routines of the EMsoft + * package, module 'dictmod' (Source/EMsoftLib/dictmod.f90) together with the + * supporting Bessel function and pseudo-random number generators of + * Source/EMsoftLib/math.f90. EMsoft is written by the Marc De Graef Research + * Group at Carnegie Mellon University. The EMsoft source attributes the original + * Expectation-Maximization implementation to Yu-Hui Chen (University of Michigan) + * and Marc De Graef (Carnegie Mellon University). + * + * Routine correspondence, EbsdLib <- EMsoft: + * + * DirectionalStats::EMforDS <- dictmod::DI_EMforDD + * DirectionalStats::Estep_ <- dictmod::DD_Estep + * DirectionalStats::Mstep_ <- dictmod::DD_Mstep + * DirectionalStats::getQandL_ <- dictmod::DD_getQandL + * DirectionalStats::Density_ <- dictmod::DD_Density + * DirectionalStats::logCp_ <- dictmod::logCp + * BesselI0 / BesselI1 / BesselIn <- math::BesselI0 / BesselI1 / BesselIn + * r8_uniform_01 <- math::r8_uniform_01 + * r8vec_uniform_01 <- math::r8vec_uniform_01 + * r8vec_normal_01 <- math::r8vec_normal_01 + * + * The 'VMF' and 'WAT' distribution selectors correspond to the EMsoft 'Dtype' + * argument. The estimator is the modified, symmetry group invariant, von + * Mises-Fisher and axial Watson distribution described in: + * + * [1] Y.H. Chen, S.U. Park, D. Wei, G. Newstadt, M.A. Jackson, J.P. Simmons, + * M. De Graef and A.O. Hero, "A Dictionary Approach to Electron Backscatter + * Diffraction Indexing", Microscopy and Microanalysis 21(3), 739-752 (2015). + * DOI: 10.1017/S1431927615000756 + * + * [2] Y.H. Chen, D. Wei, G. Newstadt, M. De Graef, J.P. Simmons and A.O. Hero, + * "Parameter Estimation in Spherical Symmetry Groups", IEEE Signal Processing + * Letters 22(8), 1152-1155 (2015). DOI: 10.1109/LSP.2014.2387206 + * + * EMsoft is available at https://github.com/EMsoft-org/EMsoft and is distributed + * under the BSD 3-Clause license reproduced below. That notice, the list of + * conditions and the disclaimer are retained here as the license requires. + * ---------------------------------------------------------------------------- + * Copyright (c) 2014-2022, Marc De Graef Research Group/Carnegie Mellon University + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * - Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * - Neither the names of Marc De Graef, Carnegie Mellon University nor the + * names of its contributors may be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + #include "DirectionalStats.hpp" #include "EbsdLib/Orientation/Quaternion.hpp" diff --git a/Source/EbsdLib/Core/DirectionalStats.hpp b/Source/EbsdLib/Core/DirectionalStats.hpp index bfb3b7c3..10eab160 100644 --- a/Source/EbsdLib/Core/DirectionalStats.hpp +++ b/Source/EbsdLib/Core/DirectionalStats.hpp @@ -1,3 +1,101 @@ +/* ============================================================================ + * Copyright (c) 2026 BlueQuartz Software, LLC + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * Redistributions in binary form must reproduce the above copyright notice, this + * list of conditions and the following disclaimer in the documentation and/or + * other materials provided with the distribution. + * + * Neither the name of BlueQuartz Software, the US Air Force, nor the names of its + * contributors may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + +/* ============================================================================ + * THIRD-PARTY ATTRIBUTION: EMsoft + * + * This file is a C++ port of the directional statistics routines of the EMsoft + * package, module 'dictmod' (Source/EMsoftLib/dictmod.f90) together with the + * supporting Bessel function and pseudo-random number generators of + * Source/EMsoftLib/math.f90. EMsoft is written by the Marc De Graef Research + * Group at Carnegie Mellon University. The EMsoft source attributes the original + * Expectation-Maximization implementation to Yu-Hui Chen (University of Michigan) + * and Marc De Graef (Carnegie Mellon University). + * + * Routine correspondence, EbsdLib <- EMsoft: + * + * DirectionalStats::EMforDS <- dictmod::DI_EMforDD + * DirectionalStats::Estep_ <- dictmod::DD_Estep + * DirectionalStats::Mstep_ <- dictmod::DD_Mstep + * DirectionalStats::getQandL_ <- dictmod::DD_getQandL + * DirectionalStats::Density_ <- dictmod::DD_Density + * DirectionalStats::logCp_ <- dictmod::logCp + * BesselI0 / BesselI1 / BesselIn <- math::BesselI0 / BesselI1 / BesselIn + * r8_uniform_01 <- math::r8_uniform_01 + * r8vec_uniform_01 <- math::r8vec_uniform_01 + * r8vec_normal_01 <- math::r8vec_normal_01 + * + * The 'VMF' and 'WAT' distribution selectors correspond to the EMsoft 'Dtype' + * argument. The estimator is the modified, symmetry group invariant, von + * Mises-Fisher and axial Watson distribution described in: + * + * [1] Y.H. Chen, S.U. Park, D. Wei, G. Newstadt, M.A. Jackson, J.P. Simmons, + * M. De Graef and A.O. Hero, "A Dictionary Approach to Electron Backscatter + * Diffraction Indexing", Microscopy and Microanalysis 21(3), 739-752 (2015). + * DOI: 10.1017/S1431927615000756 + * + * [2] Y.H. Chen, D. Wei, G. Newstadt, M. De Graef, J.P. Simmons and A.O. Hero, + * "Parameter Estimation in Spherical Symmetry Groups", IEEE Signal Processing + * Letters 22(8), 1152-1155 (2015). DOI: 10.1109/LSP.2014.2387206 + * + * EMsoft is available at https://github.com/EMsoft-org/EMsoft and is distributed + * under the BSD 3-Clause license reproduced below. That notice, the list of + * conditions and the disclaimer are retained here as the license requires. + * ---------------------------------------------------------------------------- + * Copyright (c) 2014-2022, Marc De Graef Research Group/Carnegie Mellon University + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * - Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * - Neither the names of Marc De Graef, Carnegie Mellon University nor the + * names of its contributors may be used to endorse or promote products + * derived from this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + #pragma once #include "EbsdLib/Core/EbsdLibConstants.h" diff --git a/Source/EbsdLib/Core/OrientationTransformation.hpp b/Source/EbsdLib/Core/OrientationTransformation.hpp index 220a6c1a..a07b2692 100644 --- a/Source/EbsdLib/Core/OrientationTransformation.hpp +++ b/Source/EbsdLib/Core/OrientationTransformation.hpp @@ -1159,12 +1159,19 @@ OutputType ho2ax(const InputType& h) InputType hn = h; OutputValueType sqrRtHMag = static_cast(1.0 / sqrt(hmag)); ArrayHelpers::scalarMultiply(hn, sqrRtHMag); // In place scalar multiply + if(hmag > static_cast(LPs::R1 * LPs::R1)) + { + hmag = static_cast(LPs::R1 * LPs::R1); + hm = hmag; + } + // The tfit series is valid only for 0 <= |h|^2 <= R1^2. OutputValueType s = static_cast(LPs::tfit[0] + LPs::tfit[1] * hmag); for(int i = 2; i < 16; i++) { hm = hm * hmag; s = static_cast(s + LPs::tfit[i] * hm); } + s = std::clamp(s, static_cast(-1.0), static_cast(1.0)); s = static_cast(2.0 * acos(s)); res[0] = hn[0]; res[1] = hn[1]; diff --git a/Source/EbsdLib/LaueOps/CubicLowOps.cpp b/Source/EbsdLib/LaueOps/CubicLowOps.cpp index dbd80dfc..1e7f0f17 100644 --- a/Source/EbsdLib/LaueOps/CubicLowOps.cpp +++ b/Source/EbsdLib/LaueOps/CubicLowOps.cpp @@ -87,8 +87,9 @@ constexpr std::array k_OdfNumBins = {36, 36, 36}; // Represents a 5De static const std::array k_OdfDimInitValue = {std::pow((0.75 * (ebsdlib::constants::k_PiOver2D - std::sin(ebsdlib::constants::k_PiOver2D))), (1.0 / 3.0)), std::pow((0.75 * (ebsdlib::constants::k_PiOver2D - std::sin(ebsdlib::constants::k_PiOver2D))), (1.0 / 3.0)), std::pow((0.75 * (ebsdlib::constants::k_PiOver2D - std::sin(ebsdlib::constants::k_PiOver2D))), (1.0 / 3.0))}; -static const std::array k_OdfDimStepValue = {k_OdfDimInitValue[0] / static_cast(k_OdfNumBins[0]) / 2.0, k_OdfDimInitValue[1] / static_cast(k_OdfNumBins[1]) / 2.0, - k_OdfDimInitValue[2] / static_cast(k_OdfNumBins[2]) / 2.0}; +// Dividing by half the bin count makes the inverse grid span from -init to +init, as required by the forward ODF map. +static const std::array k_OdfDimStepValue = {k_OdfDimInitValue[0] / static_cast(k_OdfNumBins[0] / 2), k_OdfDimInitValue[1] / static_cast(k_OdfNumBins[1] / 2), + k_OdfDimInitValue[2] / static_cast(k_OdfNumBins[2] / 2)}; constexpr int k_SymSize0 = 6; constexpr int k_SymSize1 = 12; @@ -191,7 +192,10 @@ constexpr double k_EtaMax = 90.0; } // namespace CubicLow // ----------------------------------------------------------------------------- -CubicLowOps::CubicLowOps() = default; +CubicLowOps::CubicLowOps() +: LaueOps(CubicLow::k_OdfDimInitValue, CubicLow::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- CubicLowOps::~CubicLowOps() = default; @@ -398,7 +402,7 @@ EulerDType CubicLowOps::determineEulerAngles(double random[3], int choose) const phi[1] = static_cast((choose / CubicLow::k_OdfNumBins[0]) % CubicLow::k_OdfNumBins[1]); phi[2] = static_cast(choose / (CubicLow::k_OdfNumBins[0] * CubicLow::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/CubicLowOps.h b/Source/EbsdLib/LaueOps/CubicLowOps.h index c81736c8..dfbad910 100644 --- a/Source/EbsdLib/LaueOps/CubicLowOps.h +++ b/Source/EbsdLib/LaueOps/CubicLowOps.h @@ -189,6 +189,7 @@ class EbsdLib_EXPORT CubicLowOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& synea) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/CubicOps.cpp b/Source/EbsdLib/LaueOps/CubicOps.cpp index c4624f1d..7a196e0e 100644 --- a/Source/EbsdLib/LaueOps/CubicOps.cpp +++ b/Source/EbsdLib/LaueOps/CubicOps.cpp @@ -278,7 +278,10 @@ constexpr double k_EtaMax = 45.0; } // namespace CubicHigh // ----------------------------------------------------------------------------- -CubicOps::CubicOps() = default; +CubicOps::CubicOps() +: LaueOps(CubicHigh::k_OdfDimInitValue, CubicHigh::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- CubicOps::~CubicOps() = default; @@ -770,7 +773,7 @@ EulerDType CubicOps::determineEulerAngles(double random[3], int choose) const phi[1] = static_cast((choose / CubicHigh::k_OdfNumBins[0]) % CubicHigh::k_OdfNumBins[1]); phi[2] = static_cast(choose / (CubicHigh::k_OdfNumBins[0] * CubicHigh::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/CubicOps.h b/Source/EbsdLib/LaueOps/CubicOps.h index c5e84c81..9aa05dc9 100644 --- a/Source/EbsdLib/LaueOps/CubicOps.h +++ b/Source/EbsdLib/LaueOps/CubicOps.h @@ -188,6 +188,7 @@ class EbsdLib_EXPORT CubicOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& synea) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp b/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp index 90f00690..1481754c 100644 --- a/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp @@ -241,7 +241,10 @@ static const SymOps k_SymOps_XParallelA = SymOps::build((choose / HexagonalLow::k_OdfNumBins[0]) % HexagonalLow::k_OdfNumBins[1]); phi[2] = static_cast(choose / (HexagonalLow::k_OdfNumBins[0] * HexagonalLow::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/HexagonalLowOps.h b/Source/EbsdLib/LaueOps/HexagonalLowOps.h index 49f70737..f684c6bf 100644 --- a/Source/EbsdLib/LaueOps/HexagonalLowOps.h +++ b/Source/EbsdLib/LaueOps/HexagonalLowOps.h @@ -189,6 +189,7 @@ class EbsdLib_EXPORT HexagonalLowOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/HexagonalOps.cpp b/Source/EbsdLib/LaueOps/HexagonalOps.cpp index a9f7458e..ae1d0fa8 100644 --- a/Source/EbsdLib/LaueOps/HexagonalOps.cpp +++ b/Source/EbsdLib/LaueOps/HexagonalOps.cpp @@ -313,7 +313,10 @@ static const SymOps k_SymOps_XParallelA = SymOps::build((choose / HexagonalHigh::k_OdfNumBins[0]) % HexagonalHigh::k_OdfNumBins[1]); phi[2] = static_cast(choose / (HexagonalHigh::k_OdfNumBins[0] * HexagonalHigh::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/HexagonalOps.h b/Source/EbsdLib/LaueOps/HexagonalOps.h index ef76dc37..9c725422 100644 --- a/Source/EbsdLib/LaueOps/HexagonalOps.h +++ b/Source/EbsdLib/LaueOps/HexagonalOps.h @@ -189,6 +189,7 @@ class EbsdLib_EXPORT HexagonalOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/LaueOps.cpp b/Source/EbsdLib/LaueOps/LaueOps.cpp index d7e649fa..c0d62837 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.cpp +++ b/Source/EbsdLib/LaueOps/LaueOps.cpp @@ -58,6 +58,9 @@ #include // for std::max #include +#include +#include +#include #include #include #include @@ -97,14 +100,110 @@ constexpr std::underlying_type_t to_underlying(Enum e) noexcept constexpr float k_OdfBinStepSize = 5.0f; +uint64_t MixBits(uint64_t value) +{ + value = (value ^ (value >> 30U)) * 0xBF58476D1CE4E5B9ULL; + value = (value ^ (value >> 27U)) * 0x94D049BB133111EBULL; + return value ^ (value >> 31U); +} + +uint64_t NextSplitMix64(uint64_t& state) +{ + state += 0x9E3779B97F4A7C15ULL; + return MixBits(state); +} + +double NextUnitInterval(uint64_t& state) +{ + return static_cast(NextSplitMix64(state) >> 11U) * 0x1.0p-53; +} + +uint64_t DoubleBits(double value) +{ + static_assert(sizeof(uint64_t) == sizeof(double)); + uint64_t bits = 0; + std::memcpy(&bits, &value, sizeof(value)); + return bits; +} + } // namespace // ----------------------------------------------------------------------------- -LaueOps::LaueOps() = default; +LaueOps::LaueOps(const std::array& odfDimInit, const std::array& odfDimStep) +: m_OdfDimInit(odfDimInit) +, m_OdfDimStep(odfDimStep) +{ +} // ----------------------------------------------------------------------------- LaueOps::~LaueOps() = default; +double LaueOps::odfBinInBallFraction(int bin) const +{ + if(bin < 0 || static_cast(bin) >= getODFSize()) + { + return 0.0; + } + + const auto bins = getOdfNumBins(); + const auto binIndex = static_cast(bin); + const std::array phi = {binIndex % bins[0], (binIndex / bins[0]) % bins[1], binIndex / (bins[0] * bins[1])}; + std::array origin; + double nearestSquared = 0.0; + double farthestSquared = 0.0; + for(size_t axis = 0; axis < 3; ++axis) + { + origin[axis] = m_OdfDimStep[axis] * static_cast(phi[axis]); + const double lower = origin[axis] - m_OdfDimInit[axis]; + const double upper = origin[axis] + m_OdfDimStep[axis] - m_OdfDimInit[axis]; + const double nearest = std::max({lower, -upper, 0.0}); + const double farthest = std::max(std::abs(lower), std::abs(upper)); + nearestSquared += nearest * nearest; + farthestSquared += farthest * farthest; + } + constexpr double k_RadiusSquared = LPs::R1 * LPs::R1; + if(nearestSquared >= k_RadiusSquared) + { + return 0.0; + } + if(farthestSquared <= k_RadiusSquared) + { + return 1.0; + } + + constexpr int k_Subdivisions = 6; + int insideCount = 0; + for(int z = 0; z < k_Subdivisions; ++z) + { + const double h3 = origin[2] + m_OdfDimStep[2] * ((z + 0.5) / k_Subdivisions) - m_OdfDimInit[2]; + for(int y = 0; y < k_Subdivisions; ++y) + { + const double h2 = origin[1] + m_OdfDimStep[1] * ((y + 0.5) / k_Subdivisions) - m_OdfDimInit[1]; + for(int x = 0; x < k_Subdivisions; ++x) + { + const double h1 = origin[0] + m_OdfDimStep[0] * ((x + 0.5) / k_Subdivisions) - m_OdfDimInit[0]; + insideCount += h1 * h1 + h2 * h2 + h3 * h3 <= k_RadiusSquared; + } + } + } + return static_cast(insideCount) / (k_Subdivisions * k_Subdivisions * k_Subdivisions); +} + +bool LaueOps::isOdfBinReachable(int bin) const +{ + return odfBinInBallFraction(bin) > 0.0; +} + +uint64_t LaueOps::clampFallbackCount() const +{ + return m_ClampFallbackCount.load(std::memory_order_relaxed); +} + +void LaueOps::resetClampFallbackCount() const +{ + m_ClampFallbackCount.store(0, std::memory_order_relaxed); +} + // ----------------------------------------------------------------------------- std::array LaueOps::getOdfBinStepSize() const { @@ -690,6 +789,44 @@ void LaueOps::_calcDetermineHomochoricValues(double random[3], double init[3], d r3 = (step[2] * phi[2]) + (step[2] * random[2]) - (init[2]); } +// ----------------------------------------------------------------------------- +bool LaueOps::_calcDetermineHomochoricValuesInBall(double random[3], double init[3], double step[3], int32_t phi[3], double& r1, double& r2, double& r3) const +{ + constexpr double k_RadiusSquared = LPs::R1 * LPs::R1; + constexpr int32_t k_MaxRedraws = 4096; + + _calcDetermineHomochoricValues(random, init, step, phi, r1, r2, r3); + auto isInBall = [&r1, &r2, &r3, k_RadiusSquared]() { return r1 * r1 + r2 * r2 + r3 * r3 <= k_RadiusSquared; }; + if(isInBall()) + { + return true; + } + + uint64_t state = 0x243F6A8885A308D3ULL; + for(size_t index = 0; index < 3; index++) + { + state = MixBits(state ^ DoubleBits(random[index])); + state = MixBits(state ^ static_cast(static_cast(phi[index]))); + } + + for(int32_t redraw = 0; redraw < k_MaxRedraws; redraw++) + { + double redrawnRandom[3] = {NextUnitInterval(state), NextUnitInterval(state), NextUnitInterval(state)}; + _calcDetermineHomochoricValues(redrawnRandom, init, step, phi, r1, r2, r3); + if(isInBall()) + { + return true; + } + } + + m_ClampFallbackCount.fetch_add(1, std::memory_order_relaxed); + const double scale = LPs::R1 / std::sqrt(r1 * r1 + r2 * r2 + r3 * r3); + r1 *= scale; + r2 *= scale; + r3 *= scale; + return false; +} + // ----------------------------------------------------------------------------- int LaueOps::_calcODFBin(double dim[3], double bins[3], double step[3], const HomochoricDType& ho) const { @@ -745,8 +882,6 @@ std::vector LaueOps::GetAllOrientationOps() /*[9]*/ m_OrientationOps.push_back(TrigonalLowOps::New()); // Trigonal-low /*[10]*/ m_OrientationOps.push_back(TrigonalOps::New()); // Trigonal-High - /*[11]*/ m_OrientationOps.push_back(OrthoRhombicOps::New()); // Axis OrthorhombicOps - return m_OrientationOps; } @@ -839,6 +974,21 @@ size_t LaueOps::getRandomSymmetryOperatorIndex(const int numSymOps) const return symOp; } +// ----------------------------------------------------------------------------- +size_t LaueOps::getRandomSymmetryOperatorIndex(const int numSymOps, std::mt19937_64& generator) const +{ + std::uniform_int_distribution distribution(0, static_cast(numSymOps - 1)); + return distribution(generator); +} + +// ----------------------------------------------------------------------------- +EulerDType LaueOps::randomizeEulerAngles(const EulerDType& euler, std::mt19937_64& generator) const +{ + const size_t symOp = getRandomSymmetryOperatorIndex(static_cast(getNumSymOps()), generator); + const QuatD qc = getQuatSymOp(symOp) * euler.toQuaternion(); + return QuaternionDType(qc).toEuler(); +} + // ----------------------------------------------------------------------------- LaueOps::Pointer LaueOps::NullPointer() { diff --git a/Source/EbsdLib/LaueOps/LaueOps.h b/Source/EbsdLib/LaueOps/LaueOps.h index fed35722..e6b39f92 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.h +++ b/Source/EbsdLib/LaueOps/LaueOps.h @@ -35,7 +35,9 @@ #pragma once #include +#include #include +#include #include #include #include @@ -110,6 +112,35 @@ class EbsdLib_EXPORT LaueOps */ virtual size_t getODFSize() const = 0; + /** + * @brief Estimates the fraction of an ODF bin inside the homochoric ball. + * @param bin Linear ODF bin index. + * @return Fraction in [0, 1], or zero for an invalid index. + * @note Boundary bins use 216 cell centers on a fixed 6-by-6-by-6 grid. + * Fully inside and fully outside boxes use exact corner bounds. A sub-grid can miss a thin intersection. + */ + double odfBinInBallFraction(int bin) const; + + /** + * @brief Tests whether the bin has a positive estimated in-ball volume. + * @param bin Linear ODF bin index. + * @return True if odfBinInBallFraction(bin) is greater than zero. + */ + bool isOdfBinReachable(int bin) const; + + /** + * @brief Returns this instance's cumulative homochoric clamp fallback count. + * @return Number of exhausted rejection-sampling attempts since the last reset. + * @note Concurrent samplers update the counter atomically. + */ + uint64_t clampFallbackCount() const; + + /** + * @brief Resets this instance's diagnostic count of homochoric clamp fallbacks. + * @note Call before sampling begins to measure a complete batch. + */ + void resetClampFallbackCount() const; + /** * @brief getNumSymmetry Returns the internal variables for k_SymSize0, k_SymSize1, k_SymSize2 * @return @@ -252,8 +283,30 @@ class EbsdLib_EXPORT LaueOps virtual EulerDType randomizeEulerAngles(const EulerDType& euler) const = 0; + /** + * @brief Selects a random symmetry operator with a clock-seeded generator. + * @param numSymOps Number of symmetry operators in the selection range. + * @return Zero-based symmetry operator index. + * @note This overload seeds a generator from the clock for each call. Use the generator-taking overload for reproducible results. + */ virtual size_t getRandomSymmetryOperatorIndex(int numSymOps) const; + /** + * @brief Selects a random symmetry operator with the specified generator. + * @param numSymOps Number of symmetry operators in the selection range. + * @param generator Generator that supplies the random stream. + * @return Zero-based symmetry operator index. + */ + size_t getRandomSymmetryOperatorIndex(int numSymOps, std::mt19937_64& generator) const; + + /** + * @brief Applies a random symmetry-equivalent rotation with the specified generator. + * @param euler Source Euler angles. + * @param generator Generator that selects the symmetry operator. + * @return Symmetry-equivalent Euler angles. + */ + EulerDType randomizeEulerAngles(const EulerDType& euler, std::mt19937_64& generator) const; + virtual RodriguesDType determineRodriguesVector(double random[3], int choose) const = 0; virtual int getOdfBin(const RodriguesDType& rod) const = 0; @@ -531,7 +584,12 @@ class EbsdLib_EXPORT LaueOps virtual bool isInsideFZ(const RodriguesDType& rod) const = 0; protected: - LaueOps(); + /** + * @brief Stores the same homochoric grid dimensions used by the derived sampler. + * @param odfDimInit Half-width of the ODF grid along each axis. + * @param odfDimStep Width of a bin along each axis. + */ + LaueOps(const std::array& odfDimInit, const std::array& odfDimStep); /** * @brief Shared annotation scaffolding for IPF images. Creates a canvas, @@ -599,6 +657,21 @@ class EbsdLib_EXPORT LaueOps */ void _calcDetermineHomochoricValues(double random[3], double init[3], double step[3], int32_t phi[3], double& r1, double& r2, double& r3) const; + /** + * @brief Samples a homochoric point in the selected ODF bin and restricts the result to the homochoric ball. + * @param random Initial offsets in the selected bin. + * @param init Half-width of each ODF grid dimension. + * @param step Width of one bin in each ODF grid dimension. + * @param phi Three-dimensional index of the selected bin. + * @param r1 Receives the first homochoric coordinate. + * @param r2 Receives the second homochoric coordinate. + * @param r3 Receives the third homochoric coordinate. + * @return True if the initial point or a redrawn point is in the homochoric ball. False if the function clamps the final point. + * @note The sampler permits 4096 redraws. The clamp is a numerical last resort that a fraction-weighted ODF should never trigger. + * Each fallback increments clampFallbackCount(). + */ + bool _calcDetermineHomochoricValuesInBall(double random[3], double init[3], double step[3], int32_t phi[3], double& r1, double& r2, double& r3) const; + /** * @brief * @param dim @@ -623,6 +696,11 @@ class EbsdLib_EXPORT LaueOps */ static QuatD ConvertToFZ(const std::vector& quatsym, const QuatD& qr, FZType fzType, AxisOrderingType order); +private: + const std::array m_OdfDimInit; + const std::array m_OdfDimStep; + mutable std::atomic m_ClampFallbackCount = 0; + public: LaueOps(const LaueOps&) = delete; // Copy Constructor Not Implemented LaueOps(LaueOps&&) = delete; // Move Constructor Not Implemented diff --git a/Source/EbsdLib/LaueOps/MonoclinicOps.cpp b/Source/EbsdLib/LaueOps/MonoclinicOps.cpp index 47746298..089798b6 100644 --- a/Source/EbsdLib/LaueOps/MonoclinicOps.cpp +++ b/Source/EbsdLib/LaueOps/MonoclinicOps.cpp @@ -87,8 +87,9 @@ constexpr std::array k_OdfNumBins = {72, 36, 72}; // Represents a 5De static const std::array k_OdfDimInitValue = {std::pow((0.7f * ((ebsdlib::constants::k_PiD)-std::sin((ebsdlib::constants::k_PiD)))), (1.0 / 3.0)), std::pow((0.75 * ((ebsdlib::constants::k_PiOver2D)-std::sin((ebsdlib::constants::k_PiOver2D)))), (1.0 / 3.0)), std::pow((0.75 * ((ebsdlib::constants::k_PiD)-std::sin((ebsdlib::constants::k_PiD)))), (1.0 / 3.0))}; -static const std::array k_OdfDimStepValue = {k_OdfDimInitValue[0] / static_cast(k_OdfNumBins[0]) / 2.0, k_OdfDimInitValue[1] / static_cast(k_OdfNumBins[1]) / 2.0, - k_OdfDimInitValue[2] / static_cast(k_OdfNumBins[2]) / 2.0}; +// Dividing by half the bin count makes the inverse grid span from -init to +init, as required by the forward ODF map. +static const std::array k_OdfDimStepValue = {k_OdfDimInitValue[0] / static_cast(k_OdfNumBins[0] / 2), k_OdfDimInitValue[1] / static_cast(k_OdfNumBins[1] / 2), + k_OdfDimInitValue[2] / static_cast(k_OdfNumBins[2] / 2)}; constexpr int k_SymSize0 = 2; constexpr int k_SymSize1 = 2; @@ -128,7 +129,10 @@ constexpr double k_ChiMax = 90.0; } // namespace Monoclinic // ----------------------------------------------------------------------------- -MonoclinicOps::MonoclinicOps() = default; +MonoclinicOps::MonoclinicOps() +: LaueOps(Monoclinic::k_OdfDimInitValue, Monoclinic::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- MonoclinicOps::~MonoclinicOps() = default; @@ -334,7 +338,7 @@ EulerDType MonoclinicOps::determineEulerAngles(double random[3], int choose) con phi[1] = static_cast((choose / Monoclinic::k_OdfNumBins[0]) % Monoclinic::k_OdfNumBins[1]); phi[2] = static_cast(choose / (Monoclinic::k_OdfNumBins[0] * Monoclinic::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/MonoclinicOps.h b/Source/EbsdLib/LaueOps/MonoclinicOps.h index 9b7ba27a..ffde9fc6 100644 --- a/Source/EbsdLib/LaueOps/MonoclinicOps.h +++ b/Source/EbsdLib/LaueOps/MonoclinicOps.h @@ -188,6 +188,7 @@ class EbsdLib_EXPORT MonoclinicOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp b/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp index b8a75122..0271e5da 100644 --- a/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp +++ b/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp @@ -140,7 +140,10 @@ constexpr double k_ChiMax = 90.0; } // namespace OrthoRhombic // ----------------------------------------------------------------------------- -OrthoRhombicOps::OrthoRhombicOps() = default; +OrthoRhombicOps::OrthoRhombicOps() +: LaueOps(OrthoRhombic::k_OdfDimInitValue, OrthoRhombic::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- OrthoRhombicOps::~OrthoRhombicOps() = default; @@ -325,7 +328,7 @@ EulerDType OrthoRhombicOps::determineEulerAngles(double random[3], int choose) c phi[1] = static_cast((choose / OrthoRhombic::k_OdfNumBins[0]) % OrthoRhombic::k_OdfNumBins[1]); phi[2] = static_cast(choose / (OrthoRhombic::k_OdfNumBins[0] * OrthoRhombic::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/OrthoRhombicOps.h b/Source/EbsdLib/LaueOps/OrthoRhombicOps.h index 75c30ea4..a017c703 100644 --- a/Source/EbsdLib/LaueOps/OrthoRhombicOps.h +++ b/Source/EbsdLib/LaueOps/OrthoRhombicOps.h @@ -190,6 +190,7 @@ class EbsdLib_EXPORT OrthoRhombicOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp b/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp index 5665aaa9..d6c772b5 100644 --- a/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp @@ -147,7 +147,10 @@ constexpr double k_ChiMax = 90.0; } // namespace TetragonalLow // ----------------------------------------------------------------------------- -TetragonalLowOps::TetragonalLowOps() = default; +TetragonalLowOps::TetragonalLowOps() +: LaueOps(TetragonalLow::k_OdfDimInitValue, TetragonalLow::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- TetragonalLowOps::~TetragonalLowOps() = default; @@ -352,7 +355,7 @@ EulerDType TetragonalLowOps::determineEulerAngles(double random[3], int choose) phi[1] = static_cast((choose / TetragonalLow::k_OdfNumBins[0]) % TetragonalLow::k_OdfNumBins[1]); phi[2] = static_cast(choose / (TetragonalLow::k_OdfNumBins[0] * TetragonalLow::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/TetragonalLowOps.h b/Source/EbsdLib/LaueOps/TetragonalLowOps.h index 7a7f3945..972bd391 100644 --- a/Source/EbsdLib/LaueOps/TetragonalLowOps.h +++ b/Source/EbsdLib/LaueOps/TetragonalLowOps.h @@ -190,6 +190,7 @@ class EbsdLib_EXPORT TetragonalLowOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/TetragonalOps.cpp b/Source/EbsdLib/LaueOps/TetragonalOps.cpp index 5a60b23e..425fc74a 100644 --- a/Source/EbsdLib/LaueOps/TetragonalOps.cpp +++ b/Source/EbsdLib/LaueOps/TetragonalOps.cpp @@ -167,7 +167,10 @@ constexpr double k_ChiMax = 90.0; } // namespace TetragonalHigh // ----------------------------------------------------------------------------- -TetragonalOps::TetragonalOps() = default; +TetragonalOps::TetragonalOps() +: LaueOps(TetragonalHigh::k_OdfDimInitValue, TetragonalHigh::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- TetragonalOps::~TetragonalOps() = default; @@ -365,7 +368,7 @@ EulerDType TetragonalOps::determineEulerAngles(double random[3], int choose) con phi[1] = static_cast((choose / TetragonalHigh::k_OdfNumBins[0]) % TetragonalHigh::k_OdfNumBins[1]); phi[2] = static_cast(choose / (TetragonalHigh::k_OdfNumBins[0] * TetragonalHigh::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/TetragonalOps.h b/Source/EbsdLib/LaueOps/TetragonalOps.h index 24e8a73a..5bb7e12e 100644 --- a/Source/EbsdLib/LaueOps/TetragonalOps.h +++ b/Source/EbsdLib/LaueOps/TetragonalOps.h @@ -190,6 +190,7 @@ class EbsdLib_EXPORT TetragonalOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/TriclinicOps.cpp b/Source/EbsdLib/LaueOps/TriclinicOps.cpp index c4b8abb6..c278b289 100644 --- a/Source/EbsdLib/LaueOps/TriclinicOps.cpp +++ b/Source/EbsdLib/LaueOps/TriclinicOps.cpp @@ -125,7 +125,10 @@ constexpr double k_ChiMax = 90.0; } // namespace Triclinic // ----------------------------------------------------------------------------- -TriclinicOps::TriclinicOps() = default; +TriclinicOps::TriclinicOps() +: LaueOps(Triclinic::k_OdfDimInitValue, Triclinic::k_OdfDimStepValue) +{ +} // ----------------------------------------------------------------------------- TriclinicOps::~TriclinicOps() = default; @@ -325,7 +328,7 @@ EulerDType TriclinicOps::determineEulerAngles(double random[3], int choose) cons phi[1] = static_cast((choose / Triclinic::k_OdfNumBins[0]) % Triclinic::k_OdfNumBins[1]); phi[2] = static_cast(choose / (Triclinic::k_OdfNumBins[0] * Triclinic::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/TriclinicOps.h b/Source/EbsdLib/LaueOps/TriclinicOps.h index 4929aff8..228b9503 100644 --- a/Source/EbsdLib/LaueOps/TriclinicOps.h +++ b/Source/EbsdLib/LaueOps/TriclinicOps.h @@ -190,6 +190,7 @@ class EbsdLib_EXPORT TriclinicOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp b/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp index 565bab74..4111e3d7 100644 --- a/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp @@ -221,7 +221,10 @@ static const SymOps k_SymOps_XParallelA = SymOps::build((choose / TrigonalLow::k_OdfNumBins[0]) % TrigonalLow::k_OdfNumBins[1]); phi[2] = static_cast(choose / (TrigonalLow::k_OdfNumBins[0] * TrigonalLow::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/TrigonalLowOps.h b/Source/EbsdLib/LaueOps/TrigonalLowOps.h index d83dec51..4a0f87a6 100644 --- a/Source/EbsdLib/LaueOps/TrigonalLowOps.h +++ b/Source/EbsdLib/LaueOps/TrigonalLowOps.h @@ -191,6 +191,7 @@ class EbsdLib_EXPORT TrigonalLowOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/LaueOps/TrigonalOps.cpp b/Source/EbsdLib/LaueOps/TrigonalOps.cpp index 497f774c..830425d4 100644 --- a/Source/EbsdLib/LaueOps/TrigonalOps.cpp +++ b/Source/EbsdLib/LaueOps/TrigonalOps.cpp @@ -241,7 +241,10 @@ static const SymOps k_SymOps_XParallelA = SymOps::build((choose / TrigonalHigh::k_OdfNumBins[0]) % TrigonalHigh::k_OdfNumBins[1]); phi[2] = static_cast(choose / (TrigonalHigh::k_OdfNumBins[0] * TrigonalHigh::k_OdfNumBins[1])); - _calcDetermineHomochoricValues(random, init, step, phi, h1, h2, h3); + _calcDetermineHomochoricValuesInBall(random, init, step, phi, h1, h2, h3); RodriguesDType ro = HomochoricDType(h1, h2, h3).toRodrigues(); ro = getODFFZRod(ro); diff --git a/Source/EbsdLib/LaueOps/TrigonalOps.h b/Source/EbsdLib/LaueOps/TrigonalOps.h index 7e327652..12989478 100644 --- a/Source/EbsdLib/LaueOps/TrigonalOps.h +++ b/Source/EbsdLib/LaueOps/TrigonalOps.h @@ -190,6 +190,7 @@ class EbsdLib_EXPORT TrigonalOps : public LaueOps int getMisoBin(const RodriguesDType& rod) const override; bool inUnitTriangle(double eta, double chi) const override; EulerDType determineEulerAngles(double random[3], int choose) const override; + using LaueOps::randomizeEulerAngles; /* Required due to C++ name hiding rules. Keeps base class 2 argument version visible */ EulerDType randomizeEulerAngles(const EulerDType& euler) const override; RodriguesDType determineRodriguesVector(double random[3], int choose) const override; int getOdfBin(const RodriguesDType& rod) const override; diff --git a/Source/EbsdLib/Orientation/Homochoric.hpp b/Source/EbsdLib/Orientation/Homochoric.hpp index f3a87146..3480ace1 100644 --- a/Source/EbsdLib/Orientation/Homochoric.hpp +++ b/Source/EbsdLib/Orientation/Homochoric.hpp @@ -8,6 +8,7 @@ #include "EbsdLib/Orientation/OrientationMatrix.hpp" #include "EbsdLib/Utilities/ModifiedLambertProjection3D.hpp" +#include #include namespace ebsdlib @@ -120,12 +121,19 @@ class Homochoric : public OrientationBase SelfType hn(*this); OutputValueType sqrRtHMag = static_cast(1.0 / sqrt(hmag)); ArrayHelpers::scalarMultiply(hn, sqrRtHMag); // In place scalar multiply + if(hmag > static_cast(LPs::R1 * LPs::R1)) + { + hmag = static_cast(LPs::R1 * LPs::R1); + hm = hmag; + } + // The tfit series is valid only for 0 <= |h|^2 <= R1^2. OutputValueType s = static_cast(LambertParametersType::tfit[0] + LambertParametersType::tfit[1] * hmag); for(int i = 2; i < 16; i++) { hm = hm * hmag; s = static_cast(s + LPs::tfit[i] * hm); } + s = std::clamp(s, static_cast(-1.0), static_cast(1.0)); s = static_cast(2.0 * acos(s)); res[0] = hn[0]; res[1] = hn[1]; diff --git a/Source/EbsdLib/Texture/StatsGen.hpp b/Source/EbsdLib/Texture/StatsGen.hpp index 6c8394cc..ab9a6986 100644 --- a/Source/EbsdLib/Texture/StatsGen.hpp +++ b/Source/EbsdLib/Texture/StatsGen.hpp @@ -485,23 +485,44 @@ class StatsGen } /** - * @brief This method will generate MDF data for a Cubic material and - * generate 1 XY scatter plots. - * @param mdf [input] This is the input MDF data which is already computed and of lenght CubicOps::k_MdfSize - * @param x [output] X Values of the Scatter plot. This memory must already be preallocated. - * @param y [outout] Y Values of the Scatter plot. This memory must already be preallocated. - * @param npoints The number of XY points for the Scatter Plot - * @param size The number of samples of the MDF to take + * @brief Generates MDF plot data with a clock-seeded generator. + * @tparam T Plot value type. + * @tparam LaueOpsType Symmetry operations for the selected crystal structure. + * @tparam ContainerType Random-access container type for the input and output arrays. + * @param mdf Precomputed MDF data. + * @param xval Receives misorientation angles in degrees. + * @param yval Receives normalized frequencies. + * @param size Number of samples to draw from the MDF. + * @return Zero if the plot generation succeeds. + * @note This overload seeds a generator from the clock for each call. Use the generator-taking overload for reproducible results. */ template static int GenMDFPlotData(ContainerType& mdf, ContainerType& xval, ContainerType& yval, int size) { - float radtodeg = 180.0f / static_cast(M_PI); - std::random_device randomDevice; // Will be used to obtain a seed for the random number engine std::mt19937_64 generator(randomDevice()); // Standard mersenne_twister_engine seeded with rd() std::mt19937_64::result_type seed = static_cast(std::chrono::steady_clock::now().time_since_epoch().count()); generator.seed(seed); + return GenMDFPlotData(mdf, xval, yval, size, generator); + } + + /** + * @brief Generates MDF plot data with the specified generator. + * @tparam T Plot value type. + * @tparam LaueOpsType Symmetry operations for the selected crystal structure. + * @tparam ContainerType Random-access container type for the input and output arrays. + * @param mdf Precomputed MDF data. + * @param xval Receives misorientation angles in degrees. + * @param yval Receives normalized frequencies. + * @param size Number of samples to draw from the MDF. + * @param generator Generator that supplies the plot sampling stream. + * @return Zero if the plot generation succeeds. + */ + template + static int GenMDFPlotData(ContainerType& mdf, ContainerType& xval, ContainerType& yval, int size, std::mt19937_64& generator) + { + float radtodeg = 180.0f / static_cast(M_PI); + std::uniform_real_distribution<> distribution(0.0, 1.0); int err = 0; diff --git a/Source/EbsdLib/Texture/Texture.hpp b/Source/EbsdLib/Texture/Texture.hpp index 4deb4b5f..b4841e19 100644 --- a/Source/EbsdLib/Texture/Texture.hpp +++ b/Source/EbsdLib/Texture/Texture.hpp @@ -35,7 +35,9 @@ #pragma once +#include #include +#include #include #include #include @@ -73,20 +75,16 @@ class Texture using ODFTableEntries = std::vector; /** - * @brief This will calculate ODF data based on an array of weights that are - * passed in and a Crystal Structure. This is templated on the container - * type, LaueOps, and type of data. Containers that adhere to the STL Vector API - * should be usable. std::vector falls into this category. The input data for the - * euler angles is in Columnar fashion instead of row major format. - * @param e1s The first euler angles - * @param e2s The second euler angles - * @param e3s The third euler angles - * @param weights Array of weights values. - * @param sigmas Array of sigma values. - * @param normalize Should the ODF data be normalized by the totalWeight value - * before returning. - * @param odf (OUT) The ODF data that is generated from this function. - * @param numEntries (OUT) The TotalWeight value that is also calculated + * @brief Calculates an ODF from weighted orientation entries and a random-texture baseline. + * @tparam T Scalar type used for weights and normalization. + * @tparam LaueOps Symmetry operations for the selected crystal structure. + * @tparam Container Random-access output container with an STL vector interface. + * @param odfTableEntries Euler angles in radians, weights, and Gaussian spread widths in bins. + * @param normalize If true, normalizes the ODF to unit total weight. + * @return ODF bin weights in the selected crystal structure's homochoric grid. + * + * The random baseline is proportional to each bin's estimated volume inside the homochoric ball. + * Weighted entries retain their Gaussian spread. Grids entirely inside the ball retain their existing arithmetic. */ template static Container CalculateODFData(const ODFTableEntries& odfTableEntries, bool normalize) @@ -200,10 +198,30 @@ class Texture else { float remainingWeight = totalWeight - totalAddWeight; - float background = remainingWeight / static_cast(ops.getODFSize()); - for(int i = 0; i < ops.getODFSize(); i++) + std::vector inBallFractions(ops.getODFSize()); + double totalFraction = 0.0; + for(size_t binIndex = 0; binIndex < inBallFractions.size(); ++binIndex) + { + inBallFractions[binIndex] = ops.odfBinInBallFraction(static_cast(binIndex)); + totalFraction += inBallFractions[binIndex]; + } + if(totalFraction == static_cast(ops.getODFSize())) { - odf[i] += background; + // Preserve the arithmetic for grids entirely inside the ball. + float background = remainingWeight / static_cast(ops.getODFSize()); + for(int i = 0; i < ops.getODFSize(); i++) + { + odf[i] += background; + } + } + else if(totalFraction > 0.0) + { + // Equal homochoric volumes have equal random-orientation probability. + const double background = static_cast(remainingWeight) / totalFraction; + for(size_t binIndex = 0; binIndex < inBallFractions.size(); ++binIndex) + { + odf[binIndex] += background * inBallFractions[binIndex]; + } } } if(normalize) @@ -219,18 +237,45 @@ class Texture } /** - * @brief CalculateMDFData Calculates MDF (Misorientation Distribution Function) data - * @param angles The angles - * @param axes The axes - * @param weights The weights - * @param odf The ODF which has been already computed and sized correctly in another function - * @param mdf [output] The MDF array to store the data which has been preallocated already - * @param numEntries The number of elemnts in teh Angles/Axes/Weights arrays which should all the be same size or at least - * the value passed here is the minium size of all the arrays. The sizes of the ODF and MDF arrays are - * determined by calling the getODFSize and getMDFSize functions of the parameterized LaueOps class. + * @brief Calculates Misorientation Distribution Function (MDF) data with a clock-seeded generator. + * @tparam T MDF value type. + * @tparam LaueOps Symmetry operations for the selected crystal structure. + * @tparam Container Random-access container type for the input and output arrays. + * @param angles Misorientation angles in radians. + * @param axes Misorientation axes with three components for each angle. + * @param weights Multiples-of-random weights. A weight of w reserves w divided by the MDF bin count of the output mass. + * @param odf Precomputed ODF data for the selected Laue class. + * @param mdf Receives the normalized MDF data. The values always sum to 1 within floating-point rounding. + * @param numEntries Number of rows to read from angles, axes, and weights. + * @note This overload seeds a generator from the clock for each call. Use the generator-taking overload for reproducible results. */ template static void CalculateMDFData(Container& angles, Container& axes, Container& weights, const Container& odf, Container& mdf, size_t numEntries) + { + std::random_device randomDevice; + std::mt19937_64 generator(randomDevice()); + std::mt19937_64::result_type seed = static_cast(std::chrono::steady_clock::now().time_since_epoch().count()); + generator.seed(seed); + CalculateMDFData(angles, axes, weights, odf, mdf, numEntries, generator); + } + + /** + * @brief Calculates Misorientation Distribution Function (MDF) data with the specified generator. + * @tparam T MDF value type. + * @tparam LaueOps Symmetry operations for the selected crystal structure. + * @tparam Container Random-access container type for the input and output arrays. + * @param angles Misorientation angles in radians. + * @param axes Misorientation axes with three components for each angle. + * @param weights Multiples-of-random weights. A weight of w reserves w divided by the MDF bin count of the output mass. + * @param odf Precomputed ODF data for the selected Laue class. + * @param mdf Receives the normalized MDF data. The values always sum to 1 within floating-point rounding. + * @param numEntries Number of rows to read from angles, axes, and weights. + * @param generator Generator that supplies the random sampling stream. + * + * If the reserved mass exceeds the 10,000-sample budget, the function scales all reserved counts proportionally. The scaled counts use the complete budget. + */ + template + static void CalculateMDFData(Container& angles, Container& axes, Container& weights, const Container& odf, Container& mdf, size_t numEntries, std::mt19937_64& generator) { LaueOps orientationOps; @@ -238,11 +283,6 @@ class Texture const int mdfsize = orientationOps.getMDFSize(); mdf.resize(orientationOps.getMDFSize()); - // Create a Random Number generator - std::random_device randomDevice; // Will be used to obtain a seed for the random number engine - std::mt19937_64 generator(randomDevice()); // Standard mersenne_twister_engine seeded with rd() - std::mt19937_64::result_type seed = static_cast(std::chrono::steady_clock::now().time_since_epoch().count()); - generator.seed(seed); std::uniform_real_distribution<> distribution(0.0, 1.0); int mbin; @@ -250,22 +290,66 @@ class Texture float totaldensity; float random1, random2, density; - for(int i = 0; i < mdfsize; i++) - { - mdf[i] = 0.0; - } - int remainingcount = 10000; - int aSize = static_cast(numEntries); + constexpr int64_t k_SampleCount = 10000; + std::vector reservedCounts(static_cast(mdfsize), 0); + int64_t totalReservedCount = 0; + const int aSize = static_cast(numEntries); for(int i = 0; i < aSize; i++) { RodriguesDType rod = AxisAngleDType(axes[3 * i], axes[3 * i + 1], axes[3 * i + 2], angles[i]).toRodrigues(); rod = orientationOps.getMDFFZRod(rod); mbin = orientationOps.getMisoBin(rod); - mdf[mbin] = static_cast(-1 * static_cast((weights[i] / static_cast(mdfsize)) * 10000.0)); - remainingcount = static_cast(remainingcount + mdf[mbin]); + const int64_t rowReservedCount = std::max(0, static_cast((weights[i] / static_cast(mdfsize)) * static_cast(k_SampleCount))); + reservedCounts[mbin] += rowReservedCount; + totalReservedCount += rowReservedCount; } + if(totalReservedCount > k_SampleCount) + { + struct ScaledRemainder + { + size_t binIndex = 0; + double remainder = 0.0; + }; + + const double scale = static_cast(k_SampleCount) / static_cast(totalReservedCount); + std::vector scaledRemainders; + int64_t scaledTotal = 0; + for(size_t binIndex = 0; binIndex < reservedCounts.size(); binIndex++) + { + if(reservedCounts[binIndex] == 0) + { + continue; + } + + const double scaledCount = static_cast(reservedCounts[binIndex]) * scale; + const int64_t integralCount = static_cast(scaledCount); + reservedCounts[binIndex] = integralCount; + scaledTotal += integralCount; + scaledRemainders.push_back({binIndex, scaledCount - static_cast(integralCount)}); + } + + std::sort(scaledRemainders.begin(), scaledRemainders.end(), [](const ScaledRemainder& lhs, const ScaledRemainder& rhs) { + if(lhs.remainder == rhs.remainder) + { + return lhs.binIndex < rhs.binIndex; + } + return lhs.remainder > rhs.remainder; + }); + for(int64_t index = 0; index < k_SampleCount - scaledTotal; index++) + { + reservedCounts[scaledRemainders[static_cast(index)].binIndex]++; + } + totalReservedCount = k_SampleCount; + } + + for(int i = 0; i < mdfsize; i++) + { + mdf[i] = static_cast(-reservedCounts[static_cast(i)]); + } + + const int remainingcount = static_cast(k_SampleCount - totalReservedCount); for(int i = 0; i < remainingcount; i++) { random1 = static_cast(distribution(generator)); @@ -287,7 +371,7 @@ class Texture choose2 = static_cast(j); } } - // This is used to create a random Homochoric vector + // The random values define a homochoric vector. std::array randx3 = {distribution(generator), distribution(generator), distribution(generator)}; EulerDType eu = orientationOps.determineEulerAngles(randx3.data(), choose1); QuatD q1 = eu.toQuaternion(); @@ -297,7 +381,7 @@ class Texture QuatD q2 = eu.toQuaternion(); RodriguesDType ro = orientationOps.calculateMisorientation(q1, q2).toRodrigues(); - ro = orientationOps.getMDFFZRod(ro); // <==== THIS IS NOT IMPELMENTED FOR ALL LAUE CLASSES + ro = orientationOps.getMDFFZRod(ro); // This operation is not implemented for all Laue classes. mbin = orientationOps.getMisoBin(ro); if(mdf[mbin] >= 0) { @@ -308,13 +392,18 @@ class Texture i = i - 1; } } + double actualTotal = 0.0; for(int i = 0; i < mdfsize; i++) { if(mdf[i] < 0) { mdf[i] = -mdf[i]; } - mdf[i] = mdf[i] / static_cast(10000.0); + actualTotal += static_cast(mdf[i]); + } + for(int i = 0; i < mdfsize; i++) + { + mdf[i] = mdf[i] / static_cast(actualTotal); } } diff --git a/Source/EbsdLib/Utilities/ODFSectionChrome.cpp b/Source/EbsdLib/Utilities/ODFSectionChrome.cpp new file mode 100644 index 00000000..9c8de595 --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionChrome.cpp @@ -0,0 +1,233 @@ +#include "EbsdLib/Utilities/ODFSectionChrome.h" + +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/Utilities/Fonts.hpp" + +#include +#include + +#include +#include + +namespace ebsdlib +{ +std::string FormatODFSectionTitle(double angleDeg) +{ + return fmt::format("\xCF\x86\xE2\x82\x82 = {:g}\xC2\xB0", angleDeg); +} + +std::string GetODFHorizontalAxisTitle() +{ + return "\xCF\x86\xE2\x82\x81"; +} + +std::string GetODFVerticalAxisTitle() +{ + return "\xCE\xA6"; +} + +double SelectODFAxisTickInterval(double maximumDeg, int32_t pixelLength) +{ + const double tenDegreePixelSpacing = static_cast(pixelLength) / (maximumDeg / 10.0); + return tenDegreePixelSpacing >= 24.0 ? 10.0 : 20.0; +} + +std::vector GenerateODFAxisTicks(double maximumDeg, int32_t pixelLength) +{ + const double intervalDeg = SelectODFAxisTickInterval(maximumDeg, pixelLength); + std::vector ticks; + for(double angleDeg = 0.0; angleDeg < maximumDeg; angleDeg += intervalDeg) + { + ticks.push_back(angleDeg); + } + if(ticks.empty() || ticks.back() != maximumDeg) + { + ticks.push_back(maximumDeg); + } + return ticks; +} + +std::vector GenerateODFAxisLabelTicks(double maximumDeg, int32_t pixelLength, float fontSize) +{ + if(fontSize <= 0.0f) + { + fontSize = std::min(std::max(10.0f, pixelLength / 24.0f) * 0.7f, std::max(8.0f, pixelLength / 32.0f) * 0.8f); + } + const auto font = fonts::GetFiraSansRegular(); + canvas_ity::canvas context(1, 1); + context.set_font(font.data(), static_cast(font.size()), fontSize); + const auto ticks = GenerateODFAxisTicks(maximumDeg, pixelLength); + std::vector labels = {ticks.front()}; + float previousRight = context.measure_text("0"); + const float endpointLeft = pixelLength - context.measure_text(fmt::format("{:g}", maximumDeg).c_str()); + constexpr float padding = 4.0f; + if(endpointLeft < previousRight + padding) + { + return labels; + } + for(size_t tickIndex = 1; tickIndex + 1 < ticks.size(); tickIndex++) + { + const double angleDeg = ticks[tickIndex]; + // Fira's numeric advance widths include the side bearings of each rendered digit. + const float width = context.measure_text(fmt::format("{:g}", angleDeg).c_str()); + const float center = static_cast(angleDeg / maximumDeg * pixelLength); + const float left = center - width / 2; + const float right = center + width / 2; + if(left >= previousRight + padding && right + padding <= endpointLeft) + { + labels.push_back(angleDeg); + previousRight = right; + } + } + if(ticks.back() != labels.back()) + { + labels.push_back(ticks.back()); + } + return labels; +} + +float ComputeODFLeftAxisGutter(float fontSize, float tickSize) +{ + // Measure ink at a fixed font size so layout checks cannot allocate an oversized canvas. + const auto font = fonts::GetFiraSansRegular(); + canvas_ity::canvas context(640, 192); + context.set_font(font.data(), static_cast(font.size()), 64.0f); + context.set_color(canvas_ity::fill_style, 0, 0, 0, 1); + context.fill_text("\xCE\xA6", 16, 96); + context.fill_text("0123456789", 128, 96); + std::vector pixels(640 * 192 * 4); + context.get_image_data(pixels.data(), 640, 192, 640 * 4, 0, 0); + int32_t titleBottom = 96; + int32_t tickTop = 96; + for(int32_t y = 0; y < 192; ++y) + { + for(int32_t x = 0; x < 640; ++x) + { + if(pixels[(y * 640 + x) * 4 + 3] != 0) + { + if(x < 128) + { + titleBottom = std::max(titleBottom, y + 1); + } + else + { + tickTop = std::min(tickTop, y); + } + } + } + } + // Rotation maps font ascent to the left of each baseline. Padding also covers pixel rounding. + return std::ceil(fontSize + (titleBottom - 96) * fontSize / 64.0f + (96 - tickTop) * tickSize / 64.0f + 10.0f); +} + +std::string GetODFColorBarTitle(ODFValueUnits sourceUnits) +{ + return sourceUnits == ODFValueUnits::MUD ? "MUD" : "MUD (from Count-Density)"; +} + +std::array GetODFSectionPanelOrigin(const ODFSectionLayoutMetrics& layout, size_t sectionIndex) +{ + return {static_cast(sectionIndex % static_cast(layout.columns)) * layout.panelSlotWidth + layout.margin + layout.leftAxisGutter, + layout.titleHeight + static_cast(sectionIndex / static_cast(layout.columns)) * layout.panelSlotHeight + layout.fontPtSize + layout.margin}; +} + +void DrawODFSectionChrome(canvas_ity::canvas& context, const ODFSectionConfiguration& config, const PreparedODFSections& sections, const ODFSectionLayoutMetrics& layout, double minimumMUD, + double maximumMUD, const std::vector& colorBar) +{ + auto boldFont = fonts::GetLatoBold(); + auto regularFont = fonts::GetLatoRegular(); + auto tickFont = fonts::GetFiraSansRegular(); + const float fontSize = layout.fontPtSize; + const float margin = layout.margin; + const float panelWidth = static_cast(layout.panelWidth); + const float panelHeight = static_cast(layout.panelHeight); + context.set_color(canvas_ity::fill_style, 0, 0, 0, 1); + context.text_baseline = canvas_ity::alphabetic; + context.set_font(boldFont.data(), static_cast(boldFont.size()), fontSize); + context.fill_text(config.title.c_str(), margin, margin + fontSize, static_cast(layout.pageWidth) - 2 * margin); + + // Fira Sans contains the Greek letters and subscripts that the embedded Lato fonts lack. + for(size_t sectionIndex = 0; sectionIndex < sections.sectionAnglesDeg.size(); sectionIndex++) + { + const auto origin = GetODFSectionPanelOrigin(layout, sectionIndex); + const float x = origin[0]; + const float y = origin[1]; + context.set_font(tickFont.data(), static_cast(tickFont.size()), fontSize); + context.text_align = canvas_ity::start; + context.fill_text(FormatODFSectionTitle(sections.sectionAnglesDeg[sectionIndex]).c_str(), x, y - margin, panelWidth); + context.fill_rectangle(x - 1, y, 1, panelHeight + 1); + context.fill_rectangle(x, y + panelHeight, panelWidth, 1); + + // PHI increases down the page, matching the prepared section row order. + const float tickSize = layout.tickFontSize; + context.set_font(tickFont.data(), static_cast(tickFont.size()), tickSize); + const auto phi1Ticks = GenerateODFAxisTicks(sections.limits.phi1MaxDeg, layout.panelWidth); + const auto phi1LabelTicks = GenerateODFAxisLabelTicks(sections.limits.phi1MaxDeg, layout.panelWidth, tickSize); + for(const double tickDeg : phi1Ticks) + { + const float fraction = static_cast(tickDeg / sections.limits.phi1MaxDeg); + const float tickX = x + panelWidth * fraction; + context.fill_rectangle(tickX - 1, y + panelHeight, 1, 3); + if(std::binary_search(phi1LabelTicks.begin(), phi1LabelTicks.end(), tickDeg)) + { + context.text_align = tickDeg == 0.0 ? canvas_ity::start : (tickDeg == sections.limits.phi1MaxDeg ? canvas_ity::rightward : canvas_ity::center); + context.fill_text(fmt::format("{:g}", tickDeg).c_str(), tickX, y + panelHeight + tickSize + 3, panelWidth / 3); + } + } + const auto phiTicks = GenerateODFAxisTicks(sections.limits.phiMaxDeg, layout.panelHeight); + const auto phiLabelTicks = GenerateODFAxisLabelTicks(sections.limits.phiMaxDeg, layout.panelHeight, tickSize); + for(const double tickDeg : phiTicks) + { + const float fraction = static_cast(tickDeg / sections.limits.phiMaxDeg); + const float tickY = y + panelHeight * fraction; + context.fill_rectangle(x - 3, tickY, 3, 1); + if(!std::binary_search(phiLabelTicks.begin(), phiLabelTicks.end(), tickDeg)) + { + continue; + } + context.save(); + context.translate(x - 4, tickY); + context.rotate(-1.5707963267948966f); + context.text_align = tickDeg == 0.0 ? canvas_ity::rightward : (tickDeg == sections.limits.phiMaxDeg ? canvas_ity::start : canvas_ity::center); + context.fill_text(fmt::format("{:g}", tickDeg).c_str(), 0, 0, panelHeight / 3); + context.restore(); + } + context.set_font(tickFont.data(), static_cast(tickFont.size()), fontSize); + context.text_align = canvas_ity::center; + context.fill_text(GetODFHorizontalAxisTitle().c_str(), x + panelWidth / 2, y + panelHeight + 2 * fontSize + margin / 2, panelWidth); + context.save(); + context.translate(std::round(x - layout.leftAxisGutter + fontSize), std::round(y + panelHeight / 2)); + context.rotate(-1.5707963267948966f); + context.fill_text(GetODFVerticalAxisTitle().c_str(), 0, 0, panelHeight); + context.restore(); + } + + const float legendX = static_cast(layout.columns) * layout.panelSlotWidth + 2 * margin; + const float legendY = layout.titleHeight + fontSize + margin; + const float barWidth = 2 * margin; + context.draw_image(colorBar.data(), 1, layout.panelHeight, 4, legendX, legendY, barWidth, panelHeight); + context.set_font(boldFont.data(), static_cast(boldFont.size()), fontSize); + context.text_align = canvas_ity::start; + context.fill_text(GetODFColorBarTitle(config.grid.units).c_str(), legendX, legendY - margin, layout.legendWidth - margin); + const float textX = legendX + barWidth + margin; + const float textWidth = static_cast(layout.pageWidth) - textX - margin; + context.set_font(regularFont.data(), static_cast(regularFont.size()), fontSize); + context.fill_text(fmt::format("{:g}", maximumMUD).c_str(), textX, legendY + fontSize / 2, textWidth); + context.fill_text(fmt::format("{:g}", minimumMUD).c_str(), legendX, legendY + panelHeight + fontSize + margin / 2, barWidth); + + const auto laueNames = LaueOps::GetLaueNames(); + const std::array labels = {fmt::format("Phase: {}", config.phaseNumber), + fmt::format("Material: {}", config.materialName), + fmt::format("Laue: {}", laueNames[config.grid.laueOpsIndex]), + fmt::format("Input: {}", config.grid.units == ODFValueUnits::MUD ? "MUD" : "CountDensity"), + "Rendered: MUD", + fmt::format("Bins: {:g}, {:g}, {:g} deg", config.grid.spacingDeg[0], config.grid.spacingDeg[1], config.grid.spacingDeg[2]), + fmt::format("Sections: {}", sections.sectionAnglesDeg.size())}; + const float infoSize = std::min(fontSize, (static_cast(layout.pageHeight) - margin - legendY - fontSize) / 8); + context.set_font(regularFont.data(), static_cast(regularFont.size()), infoSize); + for(size_t labelIndex = 0; labelIndex < labels.size(); labelIndex++) + { + context.fill_text(labels[labelIndex].c_str(), textX, legendY + fontSize + static_cast(labelIndex + 1) * infoSize, textWidth); + } +} +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/ODFSectionChrome.h b/Source/EbsdLib/Utilities/ODFSectionChrome.h new file mode 100644 index 00000000..3ac64419 --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionChrome.h @@ -0,0 +1,91 @@ +#pragma once + +#include "EbsdLib/Utilities/ODFSectionCompositor.h" + +namespace canvas_ity +{ +class canvas; +} + +namespace ebsdlib +{ +/** + * @brief Formats the φ₂ section heading with a degree symbol. + * @param angleDeg Section angle in degrees. + * @return Heading with the angle in degrees. + */ +EbsdLib_EXPORT std::string FormatODFSectionTitle(double angleDeg); + +/** + * @brief Returns the horizontal Euler-axis title as UTF-8 bytes. + * @return The φ₁ title. + */ +EbsdLib_EXPORT std::string GetODFHorizontalAxisTitle(); + +/** + * @brief Returns the vertical Euler-axis title as UTF-8 bytes. + * @return The Φ title. + */ +EbsdLib_EXPORT std::string GetODFVerticalAxisTitle(); + +/** + * @brief Selects a tick interval from the angular range and axis length. + * @param maximumDeg Positive axis maximum in degrees. + * @param pixelLength Positive axis length in pixels. + * @return Interval of 10 or 20 degrees. + */ +EbsdLib_EXPORT double SelectODFAxisTickInterval(double maximumDeg, int32_t pixelLength); + +/** + * @brief Generates axis ticks that include zero and the exact maximum. + * @param maximumDeg Positive axis maximum in degrees. + * @param pixelLength Positive axis length in pixels. + * @return Tick angles in degrees, in increasing order. + */ +EbsdLib_EXPORT std::vector GenerateODFAxisTicks(double maximumDeg, int32_t pixelLength); + +/** + * @brief Selects numeric labels with measured text widths and four pixels of padding. + * @param maximumDeg Positive axis maximum in degrees. + * @param pixelLength Positive axis length in pixels. + * @param fontSize Tick font size in pixels, or zero to derive it from the axis length. + * @return Labeled tick angles, including both endpoints when their text fits without overlap. + */ +EbsdLib_EXPORT std::vector GenerateODFAxisLabelTicks(double maximumDeg, int32_t pixelLength, float fontSize = 0.0f); + +/** + * @brief Measures the gutter needed for separate vertical title and numeric-label regions. + * @param fontSize Axis-title font size in pixels. + * @param tickSize Numeric-label font size in pixels. + * @return Gutter width in pixels, including space for tick hashes and text padding. + */ +EbsdLib_EXPORT float ComputeODFLeftAxisGutter(float fontSize, float tickSize); + +/** + * @brief Names the rendered MUD scale and identifies Count-Density sources. + * @param sourceUnits Validated source-value units. + * @return Color-bar title for the source units. + */ +EbsdLib_EXPORT std::string GetODFColorBarTitle(ODFValueUnits sourceUnits); + +/** + * @brief Returns the upper-left data-rectangle position for one section. + * @param layout Validated page metrics. + * @param sectionIndex Zero-based section index in row-major order. + * @return Horizontal and vertical pixel coordinates. + */ +EbsdLib_EXPORT std::array GetODFSectionPanelOrigin(const ODFSectionLayoutMetrics& layout, size_t sectionIndex); + +/** + * @brief Draws embedded-font titles, axes, shared color bar, and ODF metadata over the panels. + * @param context Receives the annotations on a white page with the data panels already drawn. + * @param config Validated presentation settings and source-grid metadata. + * @param sections Prepared angles and Laue limits. + * @param layout Validated page metrics. + * @param minimumMUD Applied lower scale endpoint. + * @param maximumMUD Applied upper scale endpoint. + * @param colorBar One-column RGBA image with panelHeight rows, from maximum to minimum MUD. + */ +EbsdLib_EXPORT void DrawODFSectionChrome(canvas_ity::canvas& context, const ODFSectionConfiguration& config, const PreparedODFSections& sections, const ODFSectionLayoutMetrics& layout, + double minimumMUD, double maximumMUD, const std::vector& colorBar); +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/ODFSectionCompositor.cpp b/Source/EbsdLib/Utilities/ODFSectionCompositor.cpp new file mode 100644 index 00000000..0e1463f4 --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionCompositor.cpp @@ -0,0 +1,171 @@ +#include "EbsdLib/Utilities/ODFSectionCompositor.h" +#include "EbsdLib/Utilities/ODFSectionChrome.h" + +#include +#include +#include + +#include +#include + +namespace ebsdlib +{ +namespace +{ +ODFSectionResult Failure(int32_t code, std::string message) +{ + ODFSectionResult result; + result.errorCode = code; + result.errorMessage = std::move(message); + return result; +} + +std::array MapColor(double fraction, const std::vector& controls) +{ + fraction = std::clamp(fraction, 0.0, 1.0); + size_t upper = 0; + while(upper + 4 < controls.size() && controls[upper] < fraction) + { + upper += 4; + } + const size_t lower = upper == 0 ? 0 : upper - 4; + const double weight = upper == lower ? 0.0 : std::clamp((fraction - controls[lower]) / (controls[upper] - controls[lower]), 0.0, 1.0); + std::array color = {0, 0, 0, 255}; + for(size_t componentIndex = 0; componentIndex < 3; componentIndex++) + { + color[componentIndex] = static_cast(std::lround(255.0 * ((1.0 - weight) * controls[lower + componentIndex + 1] + weight * controls[upper + componentIndex + 1]))); + } + return color; +} +} // namespace + +ODFSectionLayoutMetrics ComputeODFSectionLayout(const ODFSectionConfiguration& config) +{ + const auto limits = GetODFEulerPlotLimits(config.grid.laueOpsIndex); + if(config.sectionWidth <= 0 || config.sectionsPerRow == 0 || config.sectionsPerRow > static_cast(std::numeric_limits::max()) || config.sectionCount < 2 || + config.sectionCount > static_cast(std::numeric_limits::max()) || limits.phi1MaxDeg <= 0.0) + { + return {}; + } + ODFSectionLayoutMetrics metrics; + metrics.panelWidth = config.sectionWidth; + metrics.panelHeight = static_cast(std::lround(static_cast(config.sectionWidth) * limits.phiMaxDeg / limits.phi1MaxDeg)); + metrics.columns = static_cast(config.sectionsPerRow); + metrics.rows = static_cast((config.sectionCount + config.sectionsPerRow - 1) / config.sectionsPerRow); + metrics.fontPtSize = std::max(10.0f, static_cast(config.sectionWidth) / 24.0f); + metrics.margin = std::max(8.0f, static_cast(config.sectionWidth) / 32.0f); + metrics.tickFontSize = std::min(metrics.fontPtSize * 0.7f, metrics.margin * 0.8f); + metrics.leftAxisGutter = ComputeODFLeftAxisGutter(metrics.fontPtSize, metrics.tickFontSize); + metrics.panelSlotWidth = static_cast(metrics.panelWidth) + metrics.leftAxisGutter + 2.0f * metrics.margin; + metrics.panelSlotHeight = static_cast(metrics.panelHeight) + 3.0f * metrics.fontPtSize + 3.0f * metrics.margin; + metrics.titleHeight = metrics.fontPtSize + 2.0f * metrics.margin; + metrics.legendWidth = std::max(static_cast(config.sectionWidth) / 3.0f, 180.0f); + const float pageWidth = std::ceil(metrics.columns * metrics.panelSlotWidth + metrics.legendWidth + 2.0f * metrics.margin); + const float pageHeight = std::ceil(metrics.titleHeight + metrics.rows * metrics.panelSlotHeight + metrics.margin); + // canvas_ity stores coordinates in unsigned shorts and multiplies page dimensions as signed integers. + if(metrics.panelHeight < 1 || pageWidth >= 65535.0f || pageHeight >= 65535.0f || static_cast(pageWidth) * pageHeight > std::numeric_limits::max()) + { + return {}; + } + metrics.pageWidth = static_cast(pageWidth); + metrics.pageHeight = static_cast(pageHeight); + return metrics; +} + +ODFSectionResult ODFSectionCompositor::generateCompositeImage(const ODFSectionConfiguration& config) const +{ + if(config.sectionWidth <= 0) + { + return Failure(-7510, fmt::format("Section width must be positive; received {}.", config.sectionWidth)); + } + if(config.sectionsPerRow == 0 || config.sectionsPerRow > static_cast(std::numeric_limits::max())) + { + return Failure(-7511, fmt::format("Sections per row must be in [1, {}]; received {}.", std::numeric_limits::max(), config.sectionsPerRow)); + } + if(config.scaleMode != ODFScaleMode::Automatic && config.scaleMode != ODFScaleMode::Manual) + { + return Failure(-7512, fmt::format("Unknown ODF scale mode {}.", static_cast(config.scaleMode))); + } + if(config.scaleMode == ODFScaleMode::Manual && + (!std::isfinite(config.manualMinimumMUD) || !std::isfinite(config.manualMaximumMUD) || config.manualMinimumMUD < 0.0 || config.manualMinimumMUD >= config.manualMaximumMUD)) + { + return Failure(-7512, fmt::format("Manual MUD range must satisfy 0 <= minimum < maximum with finite endpoints; received [{}, {}].", config.manualMinimumMUD, config.manualMaximumMUD)); + } + const auto& controls = config.colorControlPoints; + if(controls.size() < 8 || controls.size() % 4 != 0) + { + return Failure(-7513, fmt::format("At least two complete [position, r, g, b] controls are required; received {} components.", controls.size())); + } + for(size_t componentIndex = 0; componentIndex < controls.size(); componentIndex++) + { + const float value = controls[componentIndex]; + if(!std::isfinite(value) || value < 0.0f || value > 1.0f || (componentIndex >= 4 && componentIndex % 4 == 0 && value <= controls[componentIndex - 4])) + { + return Failure(-7513, fmt::format("Color control component {} is invalid ({}): components must be finite in [0, 1] and positions must increase strictly.", componentIndex, value)); + } + } + const auto layout = ComputeODFSectionLayout(config); + // Reject excessive page requests before preparation allocates the section buffers. + if(layout.pageWidth == 0 && config.sectionCount >= 2 && GetODFEulerPlotLimits(config.grid.laueOpsIndex).phi1MaxDeg > 0) + { + return Failure(-7510, fmt::format("ODF layout exceeds canvas limits or has an empty panel: width {}, sections {}, columns {}.", config.sectionWidth, config.sectionCount, config.sectionsPerRow)); + } + auto prepared = PrepareODFSections(config.grid, config.sectionCount); + if(!prepared) + { + return Failure(prepared.errorCode, prepared.errorMessage); + } + for(size_t valueIndex = 0; valueIndex < prepared.sections.mudValues.size(); valueIndex++) + { + if(!std::isfinite(prepared.sections.mudValues[valueIndex])) + { + return Failure(-7512, + fmt::format("Prepared MUD value {} at section-buffer index {} is not finite; the color scale requires finite values.", prepared.sections.mudValues[valueIndex], valueIndex)); + } + } + ODFSectionResult result; + result.appliedMinimumMUD = config.scaleMode == ODFScaleMode::Automatic ? 0.0 : config.manualMinimumMUD; + result.appliedMaximumMUD = config.scaleMode == ODFScaleMode::Automatic ? prepared.sections.maximumDisplayedMUD : config.manualMaximumMUD; + if(result.appliedMaximumMUD <= result.appliedMinimumMUD) + { + return Failure(-7512, fmt::format("Automatic MUD maximum must be positive; displayed maximum is {}.", result.appliedMaximumMUD)); + } + result.width = layout.pageWidth; + result.height = layout.pageHeight; + canvas_ity::canvas context(layout.pageWidth, layout.pageHeight); + context.set_color(canvas_ity::fill_style, 1, 1, 1, 1); + context.fill_rectangle(0, 0, static_cast(layout.pageWidth), static_cast(layout.pageHeight)); + const auto& sections = prepared.sections; + std::vector panel(static_cast(layout.panelWidth) * layout.panelHeight * 4); + const double range = result.appliedMaximumMUD - result.appliedMinimumMUD; + for(size_t sectionIndex = 0; sectionIndex < config.sectionCount; sectionIndex++) + { + // Nearest-cell expansion preserves the prepared values without adding spatial interpolation. + for(int32_t row = 0; row < layout.panelHeight; row++) + { + const size_t phiIndex = static_cast(row) * sections.phiCount / static_cast(layout.panelHeight); + for(int32_t column = 0; column < layout.panelWidth; column++) + { + const size_t phi1Index = static_cast(column) * sections.phi1Count / static_cast(layout.panelWidth); + const double value = sections.mudValues[(sectionIndex * sections.phiCount + phiIndex) * sections.phi1Count + phi1Index]; + const auto color = MapColor((value - result.appliedMinimumMUD) / range, controls); + const size_t pixelOffset = (static_cast(row) * layout.panelWidth + column) * 4; + std::copy(color.begin(), color.end(), panel.begin() + pixelOffset); + } + } + const auto origin = GetODFSectionPanelOrigin(layout, sectionIndex); + context.draw_image(panel.data(), layout.panelWidth, layout.panelHeight, layout.panelWidth * 4, origin[0], origin[1], static_cast(layout.panelWidth), static_cast(layout.panelHeight)); + } + std::vector colorBar(static_cast(layout.panelHeight) * 4); + for(int32_t row = 0; row < layout.panelHeight; row++) + { + const auto color = MapColor(layout.panelHeight == 1 ? 1.0 : 1.0 - static_cast(row) / (layout.panelHeight - 1), controls); + std::copy(color.begin(), color.end(), colorBar.begin() + static_cast(row) * 4); + } + DrawODFSectionChrome(context, config, sections, layout, result.appliedMinimumMUD, result.appliedMaximumMUD, colorBar); + result.image = UInt8ArrayType::CreateArray(static_cast(layout.pageWidth) * layout.pageHeight, {4ULL}, "ODFSectionsComposite", true); + context.get_image_data(result.image->getPointer(0), layout.pageWidth, layout.pageHeight, layout.pageWidth * 4, 0, 0); + result.sectionAnglesDeg = std::move(prepared.sections.sectionAnglesDeg); + return result; +} +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/ODFSectionCompositor.h b/Source/EbsdLib/Utilities/ODFSectionCompositor.h new file mode 100644 index 00000000..47e8f1d9 --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionCompositor.h @@ -0,0 +1,113 @@ +#pragma once + +#include "EbsdLib/Utilities/ODFSectionUtilities.h" + +namespace ebsdlib +{ +/** + * @enum ODFScaleMode + * @brief Selects the shared MUD color range. + */ +enum class ODFScaleMode : uint8_t +{ + Automatic = 0, ///< Uses zero and the maximum displayed MUD value. + Manual = 1 ///< Uses the supplied MUD endpoints and clips values outside them. +}; + +/** + * @struct ODFSectionConfiguration + * @brief Supplies the borrowed grid and presentation settings for standard ODF sections. + * + * The grid values must remain valid until generateCompositeImage() returns. + */ +struct EbsdLib_EXPORT ODFSectionConfiguration +{ + ODFGridView grid; + size_t sectionCount = 6; + size_t sectionsPerRow = 3; + int32_t sectionWidth = 512; + ODFScaleMode scaleMode = ODFScaleMode::Automatic; + double manualMinimumMUD = 0.0; + double manualMaximumMUD = 1.0; + /** + * @brief At least two [position, red, green, blue] controls, with all components in [0, 1]. + * Positions must increase strictly. Values outside the control positions use the nearest endpoint color. + */ + std::vector colorControlPoints; + std::string title; + std::string materialName; + int32_t phaseNumber = 1; +}; + +/** + * @struct ODFSectionLayoutMetrics + * @brief Contains deterministic page and panel dimensions in pixels. + * + * Zero page dimensions indicate invalid or unrepresentable layout inputs. + */ +struct EbsdLib_EXPORT ODFSectionLayoutMetrics +{ + int32_t panelWidth = 0; + int32_t panelHeight = 0; + int32_t columns = 0; + int32_t rows = 0; + float fontPtSize = 0.0f; + float margin = 0.0f; + float tickFontSize = 0.0f; + float leftAxisGutter = 0.0f; + float panelSlotWidth = 0.0f; + float panelSlotHeight = 0.0f; + float titleHeight = 0.0f; + float legendWidth = 0.0f; + int32_t pageWidth = 0; + int32_t pageHeight = 0; +}; + +/** + * @struct ODFSectionResult + * @brief Contains an owned RGBA page, section metadata, or a validation error. + */ +struct EbsdLib_EXPORT ODFSectionResult +{ + UInt8ArrayType::Pointer image; + int32_t width = 0; + int32_t height = 0; + std::vector sectionAnglesDeg; + double appliedMinimumMUD = 0.0; + double appliedMaximumMUD = 0.0; + int32_t errorCode = 0; + std::string errorMessage; + /** + * @brief Reports whether composition produced an image without an error. + * @return True if the result contains an image and a nonnegative error code. + */ + explicit operator bool() const noexcept + { + return errorCode >= 0 && image != nullptr; + } +}; + +/** + * @brief Computes the page layout from the Laue limits and panel settings. + * @param config Supplies the Laue index, section count, columns, and panel width. + * @return Pixel metrics, or zero metrics if layout inputs are invalid or exceed canvas integer limits. + */ +EbsdLib_EXPORT ODFSectionLayoutMetrics ComputeODFSectionLayout(const ODFSectionConfiguration& config); + +/** + * @class ODFSectionCompositor + * @brief Prepares standard ODF sections and draws one annotated RGBA page. + */ +class EbsdLib_EXPORT ODFSectionCompositor +{ +public: + /** + * @brief Draws panels with a shared MUD color bar and embedded-font annotations. + * @param config Supplies the borrowed source grid, color controls, scale, and labels. + * @return Owned image and applied metadata, or a negative error code with a diagnostic message. + * @note Errors -7510 through -7513 identify width/layout, columns, scale, and color-control failures. + * Grid and section-count errors retain the PrepareODFSections() error codes. + */ + ODFSectionResult generateCompositeImage(const ODFSectionConfiguration& config) const; +}; +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/ODFSectionUtilities.cpp b/Source/EbsdLib/Utilities/ODFSectionUtilities.cpp new file mode 100644 index 00000000..801bc41c --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionUtilities.cpp @@ -0,0 +1,269 @@ +/* ============================================================================ + * Copyright (c) 2009-2025 BlueQuartz Software, LLC + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * Redistributions in binary form must reproduce the above copyright notice, this + * list of conditions and the following disclaimer in the documentation and/or + * other materials provided with the distribution. + * + * Neither the name of BlueQuartz Software, the US Air Force, nor the names of its + * contributors may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * The code contained herein was partially funded by the following contracts: + * United States Air Force Prime Contract FA8650-07-D-5800 + * United States Air Force Prime Contract FA8650-10-D-5210 + * United States Prime Contract Navy N00173-07-C-2068 + * + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + +#include "ODFSectionUtilities.h" + +#include "EbsdLib/Math/EbsdLibMath.h" + +#include + +#include +#include +#include +#include +#include + +namespace +{ +using namespace ebsdlib; + +constexpr int32_t k_NullValues = -7500; +constexpr int32_t k_InvalidLaueClass = -7501; +constexpr int32_t k_InvalidGrid = -7502; +constexpr int32_t k_ValueCountMismatch = -7503; +constexpr int32_t k_InvalidValue = -7504; +constexpr int32_t k_InvalidSectionCount = -7505; + +constexpr std::array k_Limits = {{ + {360.0, 90.0, 60.0}, // Hexagonal High + {360.0, 90.0, 90.0}, // Cubic High + {360.0, 180.0, 60.0}, // Hexagonal Low + {360.0, 90.0, 180.0}, // Cubic Low + {360.0, 180.0, 360.0}, // Triclinic + {360.0, 90.0, 360.0}, // Monoclinic, EbsdLib unique-axis convention + {360.0, 90.0, 180.0}, // Orthorhombic + {360.0, 180.0, 90.0}, // Tetragonal Low + {360.0, 90.0, 90.0}, // Tetragonal High + {360.0, 180.0, 120.0}, // Trigonal Low + {360.0, 90.0, 120.0} // Trigonal High +}}; + +ODFSectionPreparationResult MakeError(int32_t code, std::string message) +{ + ODFSectionPreparationResult result; + result.errorCode = code; + result.errorMessage = std::move(message); + return result; +} + +bool CoversEulerCube(const ODFGridView& grid) +{ + constexpr std::array k_FullExtentsDeg = {360.0, 180.0, 360.0}; + constexpr double k_BinCountTolerance = 1.0e-6; + for(size_t axis = 0; axis < grid.dimensions.size(); axis++) + { + const double binCount = k_FullExtentsDeg[axis] / grid.spacingDeg[axis]; + if(grid.dimensions[axis] == 0 || std::abs(binCount - static_cast(grid.dimensions[axis])) > k_BinCountTolerance) + { + return false; + } + } + return true; +} + +size_t FullIndex(size_t phi1Index, size_t phiIndex, size_t phi2Index, const std::array& dimensions) +{ + return (phi1Index * dimensions[1] + phiIndex) * dimensions[2] + phi2Index; +} + +size_t SectionIndex(size_t sectionIndex, size_t phiIndex, size_t phi1Index, size_t phi1Count, size_t phiCount) +{ + return (sectionIndex * phiCount + phiIndex) * phi1Count + phi1Index; +} + +size_t WrapIndex(int64_t index, size_t count) +{ + const auto signedCount = static_cast(count); + const int64_t remainder = index % signedCount; + return static_cast(remainder < 0 ? remainder + signedCount : remainder); +} + +size_t CountCellCentersBelow(size_t cellCount, double spacingDeg, double maximumDeg) +{ + size_t displayedCount = 0; + while(displayedCount < cellCount && (static_cast(displayedCount) + 0.5) * spacingDeg < maximumDeg) + { + displayedCount++; + } + return displayedCount; +} +} // namespace + +namespace ebsdlib +{ +ODFEulerPlotLimits GetODFEulerPlotLimits(uint32_t laueOpsIndex) +{ + if(laueOpsIndex >= k_Limits.size()) + { + return {}; + } + return k_Limits[laueOpsIndex]; +} + +std::vector GenerateODFSectionAngles(uint32_t laueOpsIndex, size_t sectionCount) +{ + if(laueOpsIndex >= k_Limits.size() || sectionCount == 0) + { + return {}; + } + + std::vector sectionAnglesDeg(sectionCount, 0.0); + const double sectionStepDeg = k_Limits[laueOpsIndex].phi2MaxDeg / static_cast(sectionCount); + for(size_t sectionIndex = 0; sectionIndex < sectionCount; sectionIndex++) + { + sectionAnglesDeg[sectionIndex] = static_cast(sectionIndex) * sectionStepDeg; + } + return sectionAnglesDeg; +} + +ODFSectionPreparationResult PrepareODFSections(const ODFGridView& grid, size_t sectionCount) +{ + if(grid.values == nullptr) + { + return MakeError(k_NullValues, "The ODF values pointer is null. Provide a valid scalar DoubleArrayType."); + } + switch(grid.units) + { + case ODFValueUnits::MUD: + case ODFValueUnits::CountDensity: + break; + default: + return MakeError(k_InvalidGrid, fmt::format("The ODF value-units code ({}) is not supported. Use MUD (0) or CountDensity (1).", static_cast(grid.units))); + } + if(grid.laueOpsIndex >= k_Limits.size()) + { + return MakeError(k_InvalidLaueClass, fmt::format("The Laue class index ({}) is not supported. Use an EbsdLib crystal-structure index from 0 through 10.", grid.laueOpsIndex)); + } + if(sectionCount < 2) + { + return MakeError(k_InvalidSectionCount, fmt::format("The section count ({}) is not valid. Use at least 2 sections.", sectionCount)); + } + if(grid.originDeg != std::array{0.0, 0.0, 0.0}) + { + return MakeError(k_InvalidGrid, + fmt::format("The ODF grid origin ({}, {}, {}) degrees is not valid. Use a zero origin for phi1, PHI, and phi2.", grid.originDeg[0], grid.originDeg[1], grid.originDeg[2])); + } + if(!std::isfinite(grid.spacingDeg[0]) || grid.spacingDeg[0] <= 0.0 || grid.spacingDeg[1] != grid.spacingDeg[0] || grid.spacingDeg[2] != grid.spacingDeg[0]) + { + return MakeError(k_InvalidGrid, fmt::format("The ODF grid spacing ({}, {}, {}) degrees is not valid. Use equal positive spacing for phi1, PHI, and phi2.", grid.spacingDeg[0], grid.spacingDeg[1], + grid.spacingDeg[2])); + } + if(!CoversEulerCube(grid)) + { + return MakeError(k_InvalidGrid, fmt::format("The ODF grid dimensions ({}, {}, {}) with spacing {} degrees do not cover the required 360 by 180 by 360 degree Euler cube.", grid.dimensions[0], + grid.dimensions[1], grid.dimensions[2], grid.spacingDeg[0])); + } + if(grid.dimensions[2] > static_cast(std::numeric_limits::max())) + { + return MakeError(k_InvalidGrid, fmt::format("The phi2 grid dimension ({}) is too large for periodic indexing.", grid.dimensions[2])); + } + + if(grid.dimensions[0] > std::numeric_limits::max() / grid.dimensions[1] || grid.dimensions[0] * grid.dimensions[1] > std::numeric_limits::max() / grid.dimensions[2]) + { + return MakeError(k_InvalidGrid, fmt::format("The ODF grid dimensions ({}, {}, {}) are too large to index safely.", grid.dimensions[0], grid.dimensions[1], grid.dimensions[2])); + } + const size_t expectedValueCount = grid.dimensions[0] * grid.dimensions[1] * grid.dimensions[2]; + if(grid.values->getNumberOfTuples() != expectedValueCount || grid.values->getNumberOfComponents() != 1) + { + return MakeError(k_ValueCountMismatch, fmt::format("The ODF array '{}' contains {} tuples with {} components, but the grid dimensions require {} scalar tuples.", grid.values->getName(), + grid.values->getNumberOfTuples(), grid.values->getNumberOfComponents(), expectedValueCount)); + } + + std::vector fullMud(expectedValueCount, 0.0); + const double stepRad = grid.spacingDeg[0] * constants::k_PiOver180D; + for(size_t phi2Index = 0; phi2Index < grid.dimensions[2]; phi2Index++) + { + for(size_t phiIndex = 0; phiIndex < grid.dimensions[1]; phiIndex++) + { + const double phiCenterRad = (static_cast(phiIndex) + 0.5) * stepRad; + for(size_t phi1Index = 0; phi1Index < grid.dimensions[0]; phi1Index++) + { + const size_t fullIndex = FullIndex(phi1Index, phiIndex, phi2Index, grid.dimensions); + const double sourceValue = grid.values->getValue(fullIndex); + if(!std::isfinite(sourceValue) || sourceValue < 0.0) + { + return MakeError(k_InvalidValue, + fmt::format("The ODF array '{}' contains the invalid value {} at linear index {}. Values must be finite and nonnegative.", grid.values->getName(), sourceValue, fullIndex)); + } + if(grid.units == ODFValueUnits::CountDensity) + { + fullMud[fullIndex] = sourceValue * 8.0 * constants::k_PiD * constants::k_PiD / (stepRad * stepRad * stepRad * std::sin(phiCenterRad)); + } + else + { + fullMud[fullIndex] = sourceValue; + } + } + } + } + + ODFSectionPreparationResult result; + auto& prepared = result.sections; + prepared.limits = k_Limits[grid.laueOpsIndex]; + prepared.sectionAnglesDeg = GenerateODFSectionAngles(grid.laueOpsIndex, sectionCount); + prepared.phi1Count = CountCellCentersBelow(grid.dimensions[0], grid.spacingDeg[0], prepared.limits.phi1MaxDeg); + prepared.phiCount = CountCellCentersBelow(grid.dimensions[1], grid.spacingDeg[0], prepared.limits.phiMaxDeg); + if(prepared.phi1Count == 0 || prepared.phiCount == 0) + { + return MakeError(k_InvalidGrid, fmt::format("The ODF grid dimensions ({}, {}, {}) with spacing ({}, {}, {}) degrees contain no usable displayed cells for Laue class index {} " + "and plotting limits (phi1={}, PHI={}, phi2={}) degrees: displayed counts phi1={}, PHI={}. Use a finer full Euler grid.", + grid.dimensions[0], grid.dimensions[1], grid.dimensions[2], grid.spacingDeg[0], grid.spacingDeg[1], grid.spacingDeg[2], grid.laueOpsIndex, + prepared.limits.phi1MaxDeg, prepared.limits.phiMaxDeg, prepared.limits.phi2MaxDeg, prepared.phi1Count, prepared.phiCount)); + } + prepared.mudValues.resize(sectionCount * prepared.phiCount * prepared.phi1Count, 0.0); + + for(size_t sectionIndex = 0; sectionIndex < sectionCount; sectionIndex++) + { + const double firstCenterDeg = 0.5 * grid.spacingDeg[2]; + const double coordinate = (prepared.sectionAnglesDeg[sectionIndex] - firstCenterDeg) / grid.spacingDeg[2]; + const double floorCoordinate = std::floor(coordinate); + const size_t lowerPhi2Index = WrapIndex(static_cast(floorCoordinate), grid.dimensions[2]); + const size_t upperPhi2Index = (lowerPhi2Index + 1) % grid.dimensions[2]; + const double fraction = coordinate - floorCoordinate; + for(size_t phiIndex = 0; phiIndex < prepared.phiCount; phiIndex++) + { + for(size_t phi1Index = 0; phi1Index < prepared.phi1Count; phi1Index++) + { + const double lowerValue = fullMud[FullIndex(phi1Index, phiIndex, lowerPhi2Index, grid.dimensions)]; + const double upperValue = fullMud[FullIndex(phi1Index, phiIndex, upperPhi2Index, grid.dimensions)]; + const double value = (1.0 - fraction) * lowerValue + fraction * upperValue; + prepared.mudValues[SectionIndex(sectionIndex, phiIndex, phi1Index, prepared.phi1Count, prepared.phiCount)] = value; + prepared.maximumDisplayedMUD = std::max(prepared.maximumDisplayedMUD, value); + } + } + } + return result; +} +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/ODFSectionUtilities.h b/Source/EbsdLib/Utilities/ODFSectionUtilities.h new file mode 100644 index 00000000..82073e19 --- /dev/null +++ b/Source/EbsdLib/Utilities/ODFSectionUtilities.h @@ -0,0 +1,162 @@ +/* ============================================================================ + * Copyright (c) 2009-2025 BlueQuartz Software, LLC + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * Redistributions in binary form must reproduce the above copyright notice, this + * list of conditions and the following disclaimer in the documentation and/or + * other materials provided with the distribution. + * + * Neither the name of BlueQuartz Software, the US Air Force, nor the names of its + * contributors may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * The code contained herein was partially funded by the following contracts: + * United States Air Force Prime Contract FA8650-07-D-5800 + * United States Air Force Prime Contract FA8650-10-D-5210 + * United States Prime Contract Navy N00173-07-C-2068 + * + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + +#pragma once + +#include "EbsdLib/Core/EbsdDataArray.hpp" +#include "EbsdLib/Core/EbsdLibConstants.h" +#include "EbsdLib/EbsdLib.h" + +#include +#include +#include +#include +#include + +namespace ebsdlib +{ +/** + * @enum ODFValueUnits + * @brief Specifies the units of an orientation distribution function grid. + */ +enum class ODFValueUnits : uint8_t +{ + MUD = 0, ///< Specifies multiples of a uniform distribution. + CountDensity = 1 ///< Specifies probability per Euler-space cell. +}; + +/** + * @struct ODFGridView + * @brief Describes a borrowed full-cube orientation distribution function grid. + * + * The caller owns the values array and must keep it valid during section preparation. + * Flat source values use phi2-fastest row-major storage: `(phi1 * nPHI + PHI) * nphi2 + phi2`. + */ +struct EbsdLib_EXPORT ODFGridView +{ + /** @brief Source scalar values in phi2-fastest row-major order. The function does not take ownership. */ + DoubleArrayType* values = nullptr; + /** @brief Cell counts listed in phi1, PHI, and phi2 order. */ + std::array dimensions = {0, 0, 0}; + /** @brief Grid origin in degrees for phi1, PHI, and phi2. */ + std::array originDeg = {0.0, 0.0, 0.0}; + /** @brief Cell spacing in degrees for phi1, PHI, and phi2. */ + std::array spacingDeg = {0.0, 0.0, 0.0}; + /** @brief Units of each source value. */ + ODFValueUnits units = ODFValueUnits::MUD; + /** @brief EbsdLib crystal-structure index that selects Euler plot limits. */ + uint32_t laueOpsIndex = CrystalStructure::UnknownCrystalStructure; +}; + +/** + * @struct ODFEulerPlotLimits + * @brief Contains maximum Euler plot angles in degrees for one Laue class. + */ +struct EbsdLib_EXPORT ODFEulerPlotLimits +{ + /** @brief Exclusive maximum phi1 cell-center angle in degrees. */ + double phi1MaxDeg = 0.0; + /** @brief Exclusive maximum PHI cell-center angle in degrees. */ + double phiMaxDeg = 0.0; + /** @brief Periodic phi2 extent in degrees. */ + double phi2MaxDeg = 0.0; +}; + +/** + * @struct PreparedODFSections + * @brief Contains interpolated and cropped orientation distribution function sections. + */ +struct EbsdLib_EXPORT PreparedODFSections +{ + /** @brief Section-major values in MUD units. PHI is the middle axis, and phi1 is fastest. */ + std::vector mudValues; + /** @brief Exact phi2 section angles in degrees. */ + std::vector sectionAnglesDeg; + /** @brief Number of displayed phi1 cell centers. */ + size_t phi1Count = 0; + /** @brief Number of displayed PHI cell centers. */ + size_t phiCount = 0; + /** @brief Maximum value in the displayed section buffer, in MUD units. */ + double maximumDisplayedMUD = 0.0; + /** @brief Euler plot limits for the selected Laue class. */ + ODFEulerPlotLimits limits; +}; + +/** + * @struct ODFSectionPreparationResult + * @brief Contains prepared sections or one validation error. + */ +struct EbsdLib_EXPORT ODFSectionPreparationResult +{ + /** @brief Prepared output. The value is empty when preparation fails. */ + PreparedODFSections sections; + /** @brief Zero for success or a negative validation error code. */ + int32_t errorCode = 0; + /** @brief Empty for success or a diagnostic message for failure. */ + std::string errorMessage; + + /** + * @brief Reports whether section preparation succeeded. + * @return True if errorCode is not negative. + */ + explicit operator bool() const noexcept + { + return errorCode >= 0; + } +}; + +/** + * @brief Returns Euler plot limits for a supported Laue class. + * @param laueOpsIndex EbsdLib crystal-structure index in the range [0, 10]. + * @return Plot limits in degrees. Unsupported indices return zero limits. + */ +EbsdLib_EXPORT ODFEulerPlotLimits GetODFEulerPlotLimits(uint32_t laueOpsIndex); + +/** + * @brief Generates exact periodic phi2 section angles for a supported Laue class. + * @param laueOpsIndex EbsdLib crystal-structure index in the range [0, 10]. + * @param sectionCount Number of equal sections across the phi2 extent. + * @return Section angles in degrees. Invalid input returns an empty vector. + */ +EbsdLib_EXPORT std::vector GenerateODFSectionAngles(uint32_t laueOpsIndex, size_t sectionCount); + +/** + * @brief Validates an Euler grid, converts its values to MUD, and interpolates periodic phi2 sections. + * @param grid Borrowed full-cube grid with scalar values and angles in degrees. + * @param sectionCount Number of equal sections. The minimum valid value is two. + * @return Prepared sections on success or a negative error code and diagnostic message on failure. + */ +EbsdLib_EXPORT ODFSectionPreparationResult PrepareODFSections(const ODFGridView& grid, size_t sectionCount); +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index 3551cacb..82ffe249 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -5,6 +5,9 @@ set(EbsdLib_${DIR_NAME}_MOC_HDRS ) set(EbsdLib_${DIR_NAME}_HDRS + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionChrome.h + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionCompositor.h + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionUtilities.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/InversePoleFigureUtilities.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureUtilities.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ModifiedLambertProjection.h @@ -41,6 +44,9 @@ set(EbsdLib_${DIR_NAME}_HDRS ) set(EbsdLib_${DIR_NAME}_SRCS + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionChrome.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionCompositor.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ODFSectionUtilities.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/InversePoleFigureUtilities.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureUtilities.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ModifiedLambertProjection.cpp diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index b5b5ef4f..81d151b9 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -13,6 +13,9 @@ include(Catch) # ------------------------------------------------------------------------------ # Define the list of unit test source files set(EbsdLib_UnitTest_SRCS + ${EbsdLibProj_SOURCE_DIR}/Source/Test/CanvasItyImplementation.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/ODFSectionCompositorTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/ODFSectionUtilitiesTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/AngImportTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/AngleFileLoaderTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/OrientationTest.cpp diff --git a/Source/Test/CanvasItyImplementation.cpp b/Source/Test/CanvasItyImplementation.cpp new file mode 100644 index 00000000..abd0784b --- /dev/null +++ b/Source/Test/CanvasItyImplementation.cpp @@ -0,0 +1,2 @@ +#define CANVAS_ITY_IMPLEMENTATION +#include diff --git a/Source/Test/ConvertToFundamentalZoneTest.cpp b/Source/Test/ConvertToFundamentalZoneTest.cpp index 3823e4fe..2fb86496 100644 --- a/Source/Test/ConvertToFundamentalZoneTest.cpp +++ b/Source/Test/ConvertToFundamentalZoneTest.cpp @@ -82,8 +82,9 @@ static RodriguesDType convertRodrigues(const std::array& rod) TEST_CASE("ebsdlib::ConvertToFundamentalZoneTest", "[EbsdLib][ConvertToFundamentalZoneTest]") { auto ops = LaueOps::GetAllOrientationOps(); + REQUIRE(ops.size() == detail::k_FZValues.size()); std::cout << "############################################################\n"; - for(size_t opsIdx = 0; opsIdx < ops.size() - 1; ++opsIdx) // We ONLY want Cubic 432 rotation group + for(size_t opsIdx = 0; opsIdx < ops.size(); ++opsIdx) { std::cout << "OpsIndex: " << opsIdx << " " << ops[opsIdx]->getRotationPointGroup() << ", " << ops[opsIdx]->getSymmetryName() << ", " << ops[opsIdx]->getPointGroup() << ", " << ops[opsIdx]->FZTypeToString(ops[opsIdx]->getFZType()) << ", " << ops[opsIdx]->AxisOrderingTypeToString(ops[opsIdx]->getAxisOrderingType()) << std::endl; diff --git a/Source/Test/DirectionalStatsTest.cpp b/Source/Test/DirectionalStatsTest.cpp index a1f886ba..4f900a29 100644 --- a/Source/Test/DirectionalStatsTest.cpp +++ b/Source/Test/DirectionalStatsTest.cpp @@ -140,7 +140,7 @@ TEST_CASE("DirectionalStatsTest:AverageOrientation", "[DirectionalStatsTest]") for(const auto& op : ops) { const std::string rpg = op->getRotationPointGroup(); - // Skip Triclinic (no FZ boundary) and duplicates (OrthoRhombicOps appears twice) + // Skip Triclinic (no FZ boundary) and rotation point groups already tested. if(rpg == "1" || tested.count(rpg) > 0) { continue; diff --git a/Source/Test/LaueOpsTest.cpp b/Source/Test/LaueOpsTest.cpp index 5c3d1898..8ea9edc4 100644 --- a/Source/Test/LaueOpsTest.cpp +++ b/Source/Test/LaueOpsTest.cpp @@ -1,5 +1,6 @@ #include +#include "EbsdLib/Core/EbsdLibConstants.h" #include "EbsdLib/LaueOps/CubicLowOps.h" #include "EbsdLib/LaueOps/CubicOps.h" #include "EbsdLib/LaueOps/HexagonalLowOps.h" @@ -12,17 +13,38 @@ #include "EbsdLib/LaueOps/TriclinicOps.h" #include "EbsdLib/LaueOps/TrigonalLowOps.h" #include "EbsdLib/LaueOps/TrigonalOps.h" +#include "EbsdLib/Orientation/AxisAngle.hpp" +#include "EbsdLib/Orientation/Homochoric.hpp" #include "EbsdLib/Orientation/OrientationMatrix.hpp" #include "EbsdLib/Orientation/Rodrigues.hpp" #include "EbsdLib/Utilities/ColorTable.h" +#include #include +#include #include +#include +#include #include #include using namespace ebsdlib; +template +concept HasGeneratorRandomizeEulerAngles = requires(OpsType ops, const EulerDType& euler, std::mt19937_64& generator) { ops.randomizeEulerAngles(euler, generator); }; + +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); +static_assert(HasGeneratorRandomizeEulerAngles); + // ----------------------------------------------------------------------------- // getDefaultPoleFigureNames returns the plotted plane-normal families in brace // (plane-family) notation. The family identity does NOT depend on HexConvention: @@ -191,8 +213,7 @@ TEST_CASE("ebsdlib::LaueOpsTest::GenerateIPFTriangleLegend_HexConvention_Hexagon TEST_CASE("ebsdlib::LaueOpsTest::GetAllOrientationOps", "[EbsdLib][LaueOpsTest]") { auto ops = LaueOps::GetAllOrientationOps(); - // Should return exactly 12 entries (one for each Laue group index 0-11) - REQUIRE(ops.size() == 12); + REQUIRE(ops.size() == CrystalStructure::LaueGroupEnd); for(size_t i = 0; i < ops.size(); i++) { @@ -200,6 +221,94 @@ TEST_CASE("ebsdlib::LaueOpsTest::GetAllOrientationOps", "[EbsdLib][LaueOpsTest]" } } +// ----------------------------------------------------------------------------- +TEST_CASE("ebsdlib::LaueOpsTest::DetermineEulerAnglesSamplesValidOrientations", "[EbsdLib][LaueOpsTest]") +{ + constexpr double k_MaxRotationAngle = ebsdlib::constants::k_PiD + 1.0e-9; + constexpr double k_MinimumBinAgreement = 0.15; + constexpr size_t k_TargetRoundTripSamples = 4000; + uint64_t randomState = 0x4D595DF4D0F33173ULL; + + auto nextRandom = [&randomState]() { + randomState += 0x9E3779B97F4A7C15ULL; + uint64_t value = randomState; + value = (value ^ (value >> 30U)) * 0xBF58476D1CE4E5B9ULL; + value = (value ^ (value >> 27U)) * 0x94D049BB133111EBULL; + value ^= value >> 31U; + return static_cast(value >> 11U) * 0x1.0p-53; + }; + + const auto allOps = LaueOps::GetAllOrientationOps(); + for(const auto& ops : allOps) + { + const size_t roundTripStride = std::max(1, (ops->getODFSize() + k_TargetRoundTripSamples - 1) / k_TargetRoundTripSamples); + size_t matchingBinCount = 0; + size_t roundTripSampleCount = 0; + for(size_t bin = 0; bin < ops->getODFSize(); bin++) + { + double random[3] = {nextRandom(), nextRandom(), nextRandom()}; + const EulerDType first = ops->determineEulerAngles(random, static_cast(bin)); + const EulerDType second = ops->determineEulerAngles(random, static_cast(bin)); + + for(size_t component = 0; component < 3; component++) + { + if(!std::isfinite(first[component])) + { + FAIL("Non-finite Euler component for " << ops->getNameOfClass() << " at ODF bin " << bin << ", component " << component); + } + if(first[component] != second[component]) + { + FAIL("Non-deterministic Euler component for " << ops->getNameOfClass() << " at ODF bin " << bin << ", component " << component); + } + } + + const AxisAngleDType axisAngle = first.toAxisAngle(); + if(!std::isfinite(axisAngle[3]) || axisAngle[3] > k_MaxRotationAngle) + { + FAIL("Invalid rotation angle for " << ops->getNameOfClass() << " at ODF bin " << bin << ": " << axisAngle[3]); + } + + if(bin % roundTripStride == 0) + { + const int roundTripBin = ops->getOdfBin(first.toRodrigues()); + matchingBinCount += roundTripBin == static_cast(bin) ? 1 : 0; + roundTripSampleCount++; + } + } + + const double binAgreement = static_cast(matchingBinCount) / static_cast(roundTripSampleCount); + // The threshold separates a consistent grid from an inconsistent grid. It is not a measure of folding accuracy. + INFO(ops->getNameOfClass() << " sampled ODF bin agreement: " << binAgreement); + CHECK(binAgreement >= k_MinimumBinAgreement); + } +} + +// ----------------------------------------------------------------------------- +TEST_CASE("ebsdlib::LaueOpsTest::HomochoricToAxisAngleClampsOutsideDomain", "[EbsdLib][LaueOpsTest]") +{ + constexpr size_t k_NumSteps = 4096; + constexpr double k_MaxRotationAngle = ebsdlib::constants::k_PiD + 1.0e-9; + + for(size_t step = 0; step <= k_NumSteps; step++) + { + const double magnitude = 2.5 * LPs::R1 * static_cast(step) / static_cast(k_NumSteps); + const HomochoricDType homochoric(LPs::isrt * magnitude, LPs::isrt * magnitude, LPs::isrt * magnitude); + const AxisAngleDType axisAngle = homochoric.toAxisAngle(); + + for(size_t component = 0; component < 4; component++) + { + if(!std::isfinite(axisAngle[component])) + { + FAIL("Non-finite axis-angle component at homochoric magnitude " << magnitude << ", component " << component); + } + } + if(axisAngle[3] > k_MaxRotationAngle) + { + FAIL("Rotation angle exceeds pi at homochoric magnitude " << magnitude << ": " << axisAngle[3]); + } + } +} + // ----------------------------------------------------------------------------- TEST_CASE("ebsdlib::LaueOpsTest::GetNumSymOps", "[EbsdLib][LaueOpsTest]") { @@ -239,7 +348,7 @@ TEST_CASE("ebsdlib::LaueOpsTest::GetSymmetryName", "[EbsdLib][LaueOpsTest]") TEST_CASE("ebsdlib::LaueOpsTest::GetLaueNames", "[EbsdLib][LaueOpsTest]") { auto names = LaueOps::GetLaueNames(); - REQUIRE(names.size() == 12); + REQUIRE(names.size() == CrystalStructure::LaueGroupEnd); for(const auto& name : names) { @@ -547,6 +656,56 @@ bool quatsSameRotation(const QuatD& a, const QuatD& b, double tol) } } // namespace +// ----------------------------------------------------------------------------- +TEST_CASE("ebsdlib::LaueOpsTest::SeededRandomSymmetry", "[EbsdLib][LaueOpsTest]") +{ + const auto allOps = LaueOps::GetAllOrientationOps(); + const EulerDType inputEuler(0.37, 0.91, 1.42); + const QuatD inputQuat = inputEuler.toQuaternion(); + constexpr uint64_t k_Seed = 0x5EED1234ULL; + constexpr double k_QuaternionTolerance = 1.0e-12; + + for(const auto& ops : allOps) + { + INFO(ops->getNameOfClass()); + std::mt19937_64 firstGenerator(k_Seed); + std::mt19937_64 secondGenerator(k_Seed); + + for(size_t draw = 0; draw < 100; draw++) + { + const EulerDType firstEuler = ops->randomizeEulerAngles(inputEuler, firstGenerator); + const EulerDType secondEuler = ops->randomizeEulerAngles(inputEuler, secondGenerator); + for(size_t component = 0; component < 3; component++) + { + CHECK(firstEuler[component] == secondEuler[component]); + } + + const QuatD randomizedQuat = firstEuler.toQuaternion(); + bool foundEquivalentOperator = false; + for(size_t symOp = 0; symOp < ops->getNumSymOps(); symOp++) + { + const QuatD expectedQuat = ops->getQuatSymOp(symOp) * inputQuat; + if(quatsSameRotation(randomizedQuat, expectedQuat, k_QuaternionTolerance)) + { + foundEquivalentOperator = true; + break; + } + } + CHECK(foundEquivalentOperator); + } + + std::vector operatorWasSelected(ops->getNumSymOps(), false); + std::mt19937_64 indexGenerator(k_Seed); + for(size_t draw = 0; draw < 10000; draw++) + { + const size_t symOp = ops->getRandomSymmetryOperatorIndex(static_cast(ops->getNumSymOps()), indexGenerator); + REQUIRE(symOp < operatorWasSelected.size()); + operatorWasSelected[symOp] = true; + } + CHECK(std::all_of(operatorWasSelected.cbegin(), operatorWasSelected.cend(), [](bool selected) { return selected; })); + } +} + // ----------------------------------------------------------------------------- // Validates that all three symmetry operator representations (quaternion, // Rodrigues, rotation matrix) are mutually consistent at each index for diff --git a/Source/Test/ODFSectionCompositorTest.cpp b/Source/Test/ODFSectionCompositorTest.cpp new file mode 100644 index 00000000..39f821d5 --- /dev/null +++ b/Source/Test/ODFSectionCompositorTest.cpp @@ -0,0 +1,581 @@ +#include + +#include "EbsdLib/Utilities/Fonts.hpp" +#include "EbsdLib/Utilities/ODFSectionChrome.h" +#include "EbsdLib/Utilities/ODFSectionCompositor.h" + +#include + +#include +#include +#include +#include +#include + +using namespace ebsdlib; + +namespace +{ +ODFSectionConfiguration RenderConfiguration(DoubleArrayType* values) +{ + ODFSectionConfiguration config; + config.grid = {values, {4, 2, 4}, {0, 0, 0}, {90, 90, 90}, ODFValueUnits::MUD, CrystalStructure::Triclinic}; + config.sectionCount = 4; + config.sectionsPerRow = 2; + config.sectionWidth = 128; + config.scaleMode = ODFScaleMode::Manual; + config.manualMaximumMUD = 4; + config.colorControlPoints = {0, 0, 0, 1, 1, 1, 0, 0}; + config.title = "Standard ODF sections"; + config.materialName = "Titanium"; + return config; +} + +void RequirePixel(const ODFSectionResult& result, int32_t x, int32_t y, std::array expected) +{ + const size_t offset = (static_cast(y) * result.width + x) * 4; + for(size_t componentIndex = 0; componentIndex < 4; componentIndex++) + { + INFO("Pixel " << x << ", " << y << " component " << componentIndex); + REQUIRE(std::abs(static_cast(result.image->getValue(offset + componentIndex)) - expected[componentIndex]) <= 1); + } +} + +bool HasInk(const ODFSectionResult& result, int32_t x, int32_t y, int32_t width, int32_t height) +{ + for(int32_t row = y; row < y + height; row++) + { + for(int32_t column = x; column < x + width; column++) + { + if(result.image->getValue((static_cast(row) * result.width + column) * 4) < 128) + { + return true; + } + } + } + return false; +} + +bool HasPixelDifference(const ODFSectionResult& first, const ODFSectionResult& second, int32_t x, int32_t y, int32_t width, int32_t height) +{ + for(int32_t row = y; row < y + height; row++) + { + for(int32_t column = x; column < x + width; column++) + { + const size_t pixelOffset = (static_cast(row) * first.width + column) * 4; + for(size_t componentIndex = 0; componentIndex < 4; componentIndex++) + { + if(first.image->getValue(pixelOffset + componentIndex) != second.image->getValue(pixelOffset + componentIndex)) + { + return true; + } + } + } + } + return false; +} + +std::vector RenderFiraGlyph(const std::vector& font, const std::string& glyph) +{ + canvas_ity::canvas context(48, 48); + context.set_color(canvas_ity::fill_style, 1, 1, 1, 1); + context.fill_rectangle(0, 0, 48, 48); + context.set_color(canvas_ity::fill_style, 0, 0, 0, 1); + context.set_font(font.data(), static_cast(font.size()), 24); + context.fill_text(glyph.c_str(), 4, 32); + std::vector pixels(48 * 48 * 4); + context.get_image_data(pixels.data(), 48, 48, 48 * 4, 0, 0); + return pixels; +} +} // namespace + +TEST_CASE("ebsdlib::ODFSectionCompositor::ChromeTextAndTicks", "[EbsdLib][ODFSectionCompositor]") +{ + const std::string title = FormatODFSectionTitle(10.0); + REQUIRE(std::vector(title.begin(), title.end()) == std::vector{0xCF, 0x86, 0xE2, 0x82, 0x82, 0x20, 0x3D, 0x20, 0x31, 0x30, 0xC2, 0xB0}); + const std::string horizontalAxisTitle = GetODFHorizontalAxisTitle(); + REQUIRE(std::vector(horizontalAxisTitle.begin(), horizontalAxisTitle.end()) == std::vector{0xCF, 0x86, 0xE2, 0x82, 0x81}); + const std::string verticalAxisTitle = GetODFVerticalAxisTitle(); + REQUIRE(std::vector(verticalAxisTitle.begin(), verticalAxisTitle.end()) == std::vector{0xCE, 0xA6}); + REQUIRE(SelectODFAxisTickInterval(360.0, 900) == Approx(10.0)); + REQUIRE(SelectODFAxisTickInterval(360.0, 864) == Approx(10.0)); + REQUIRE(SelectODFAxisTickInterval(360.0, 512) == Approx(20.0)); + REQUIRE(GenerateODFAxisTicks(90.0, 128) == std::vector{0.0, 20.0, 40.0, 60.0, 80.0, 90.0}); + REQUIRE(GenerateODFAxisTicks(360.0, 128) == std::vector{0.0, 20.0, 40.0, 60.0, 80.0, 100.0, 120.0, 140.0, 160.0, 180.0, 200.0, 220.0, 240.0, 260.0, 280.0, 300.0, 320.0, 340.0, 360.0}); + REQUIRE(GenerateODFAxisLabelTicks(360.0, 128).size() < GenerateODFAxisTicks(360.0, 128).size()); + REQUIRE(GenerateODFAxisTicks(360.0, 864).size() == 37); + REQUIRE(GetODFColorBarTitle(ODFValueUnits::MUD) == "MUD"); + REQUIRE(GetODFColorBarTitle(ODFValueUnits::CountDensity) == "MUD (from Count-Density)"); + for(const double tick : GenerateODFAxisTicks(90.0, 128)) + { + const std::string label = fmt::format("{:g}", tick); + REQUIRE(label.find("deg") == std::string::npos); + for(const unsigned char value : label) + { + REQUIRE(value < 0x80); + } + } +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::FiraGreekGlyphs", "[EbsdLib][ODFSectionCompositor]") +{ + const auto font = fonts::GetFiraSansRegular(); + const auto missingGlyph = RenderFiraGlyph(font, "\xF4\x8F\xBF\xBF"); + const std::array requiredGlyphs = {"\xCF\x86", "\xE2\x82\x81", "\xE2\x82\x82", "\xCE\xA6", "\xC2\xB0"}; + for(const auto& glyph : requiredGlyphs) + { + const auto pixels = RenderFiraGlyph(font, glyph); + REQUIRE(pixels != missingGlyph); + REQUIRE(std::any_of(pixels.begin(), pixels.end(), [](uint8_t value) { return value < 128; })); + } +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::AxisLabelBounds", "[EbsdLib][ODFSectionCompositor]") +{ + const int32_t width = GENERATE(128, 512, 864); + const bool vertical = GENERATE(false, true); + const double maximum = vertical ? 90.0 : 360.0; + const int32_t length = vertical ? width / 4 : width; + const float tickSize = std::min(std::max(10.0f, width / 24.0f) * 0.7f, std::max(8.0f, width / 32.0f) * 0.8f); + const auto font = fonts::GetFiraSansRegular(); + canvas_ity::canvas context(1, 1); + context.set_font(font.data(), static_cast(font.size()), tickSize); + const auto labels = GenerateODFAxisLabelTicks(maximum, length, tickSize); + REQUIRE(labels.front() == 0.0); + REQUIRE(labels.back() == maximum); + float previousRight = -4.0f; + for(const double tick : labels) + { + const float textWidth = context.measure_text(fmt::format("{:g}", tick).c_str()); + const float position = static_cast(tick / maximum * length); + const float left = position - (tick == 0 ? 0 : (tick == maximum ? textWidth : textWidth / 2)); + CAPTURE(width, vertical, tick, left, previousRight); + CHECK(left >= previousRight + 4.0f); + CHECK(left >= 0.0f); + CHECK(left + textWidth <= length); + previousRight = left + textWidth; + } +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::ComposedVerticalGlyphBounds", "[EbsdLib][ODFSectionCompositor]") +{ + const int32_t width = GENERATE(128, 512, 864); + auto values = DoubleArrayType::CreateArray(32, "Uniform", true); + values->initializeWithValue(1.0); + auto config = RenderConfiguration(values.get()); + config.sectionWidth = width; + config.sectionCount = 6; + config.sectionsPerRow = 3; + config.grid.laueOpsIndex = GENERATE(CrystalStructure::Triclinic, CrystalStructure::Hexagonal_High); + const auto layout = ComputeODFSectionLayout(config); + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE(result); + const auto font = fonts::GetFiraSansRegular(); + // An isolated glyph with a white border detects clipping, missing pixels, and adjacent tick-label ink. + for(size_t sectionIndex = 0; sectionIndex < config.sectionCount; ++sectionIndex) + { + const auto origin = GetODFSectionPanelOrigin(layout, sectionIndex); + const int32_t baselineX = static_cast(std::lround((sectionIndex % config.sectionsPerRow) * layout.panelSlotWidth + layout.margin + layout.fontPtSize)); + const int32_t baselineY = static_cast(std::lround(origin[1] + layout.panelHeight / 2.0f)); + const int32_t extent = static_cast(std::ceil(layout.fontPtSize)) + 3; + canvas_ity::canvas reference(result.width, result.height); + reference.set_color(canvas_ity::fill_style, 1, 1, 1, 1); + reference.fill_rectangle(0, 0, static_cast(result.width), static_cast(result.height)); + reference.set_color(canvas_ity::fill_style, 0, 0, 0, 1); + reference.set_font(font.data(), static_cast(font.size()), layout.fontPtSize); + reference.text_align = canvas_ity::center; + reference.translate(static_cast(baselineX), static_cast(baselineY)); + reference.rotate(-1.5707963267948966f); + reference.fill_text("\xCE\xA6", 0, 0); + std::vector pixels(static_cast(4 * extent * extent * 4)); + reference.get_image_data(pixels.data(), 2 * extent, 2 * extent, 2 * extent * 4, baselineX - extent, baselineY - extent); + int32_t minX = 2 * extent, minY = 2 * extent, maxX = 0, maxY = 0; + for(int32_t y = 0; y < 2 * extent; ++y) + { + for(int32_t x = 0; x < 2 * extent; ++x) + { + if(pixels[(y * 2 * extent + x) * 4] < 250) + { + minX = std::min(minX, x); + maxX = std::max(maxX, x); + minY = std::min(minY, y); + maxY = std::max(maxY, y); + } + } + } + REQUIRE(minX < maxX); + CHECK(baselineX - extent + minX > 0); + CHECK(baselineX - extent + maxX + 2 < origin[0] - 3); + bool glyphMatches = true; + for(int32_t y = minY - 2; y <= maxY + 2; ++y) + { + for(int32_t x = minX - 2; x <= maxX + 2; ++x) + { + const int32_t pageX = baselineX - extent + x; + const int32_t pageY = baselineY - extent + y; + REQUIRE(pageX >= 0); + REQUIRE(pageX < result.width); + REQUIRE(pageY >= 0); + REQUIRE(pageY < result.height); + const int expected = pixels[(y * 2 * extent + x) * 4]; + const int actual = result.image->getValue((static_cast(pageY) * result.width + pageX) * 4); + glyphMatches = glyphMatches && std::abs(expected - actual) <= 1; + } + } + CAPTURE(width, sectionIndex); + CHECK(glyphMatches); + const auto limits = GetODFEulerPlotLimits(config.grid.laueOpsIndex); + for(const double tick : GenerateODFAxisTicks(limits.phi1MaxDeg, layout.panelWidth)) + { + const auto x = static_cast(std::ceil(origin[0] + tick / limits.phi1MaxDeg * layout.panelWidth)) - 1; + const auto y = static_cast(std::ceil(origin[1] + layout.panelHeight)) + 1; + // Fractional tick positions can cover less than half of the sampled pixel. + CHECK(result.image->getValue((static_cast(y) * result.width + x) * 4) < 255); + } + for(const double tick : GenerateODFAxisTicks(limits.phiMaxDeg, layout.panelHeight)) + { + const auto x = static_cast(std::ceil(origin[0])) - 2; + const auto y = static_cast(std::floor(origin[1] + tick / limits.phiMaxDeg * layout.panelHeight)); + CHECK(result.image->getValue((static_cast(y) * result.width + x) * 4) < 255); + } + } +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::DefaultsAndLayout", "[EbsdLib][ODFSectionCompositor]") +{ + ODFSectionConfiguration config; + REQUIRE(config.sectionCount == 6); + REQUIRE(config.sectionsPerRow == 3); + REQUIRE(config.sectionWidth == 512); + REQUIRE(config.scaleMode == ODFScaleMode::Automatic); + REQUIRE(config.manualMinimumMUD == 0.0); + REQUIRE(config.manualMaximumMUD == 1.0); + REQUIRE(config.phaseNumber == 1); + REQUIRE(config.colorControlPoints.empty()); + config.grid.laueOpsIndex = CrystalStructure::Hexagonal_High; + const auto layout = ComputeODFSectionLayout(config); + REQUIRE(layout.panelWidth == 512); + REQUIRE(layout.panelHeight == 128); + REQUIRE(layout.columns == 3); + REQUIRE(layout.rows == 2); + REQUIRE(layout.fontPtSize == Approx(512.0 / 24.0)); + REQUIRE(layout.margin == 16.0f); + REQUIRE(layout.leftAxisGutter == 41.0f); + REQUIRE(layout.panelSlotWidth == 585.0f); + REQUIRE(layout.panelSlotHeight == 240.0f); + REQUIRE(layout.titleHeight == Approx(53.333333)); + REQUIRE(layout.legendWidth == 180.0f); + REQUIRE(layout.pageWidth == 1967); + REQUIRE(layout.pageHeight == 550); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::EmptyDisplayedCrop", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(4, "CoarseGrid", true); + values->initializeWithValue(1.0); + auto config = RenderConfiguration(values.get()); + config.grid = {values.get(), {2, 1, 2}, {0, 0, 0}, {180, 180, 180}, ODFValueUnits::MUD, CrystalStructure::Hexagonal_High}; + config.sectionCount = 6; + config.sectionsPerRow = 3; + config.manualMaximumMUD = 1.0; + config.scaleMode = GENERATE(ODFScaleMode::Manual, ODFScaleMode::Automatic); + // Preparation must reject the empty crop before the renderer can read its values. + REQUIRE_FALSE(PrepareODFSections(config.grid, config.sectionCount)); + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + REQUIRE(result.image == nullptr); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::InvalidConfiguration", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(32, "Uniform", true); + values->initializeWithValue(1.0); + ODFSectionConfiguration config; + config.grid = {values.get(), {4, 2, 4}, {0, 0, 0}, {90, 90, 90}, ODFValueUnits::MUD, CrystalStructure::Triclinic}; + config.colorControlPoints = {0, 0, 0, 1, 1, 1, 0, 0}; + int32_t expectedCode = 0; + SECTION("Nonpositive width") + { + config.sectionWidth = 0; + expectedCode = -7510; + } + SECTION("Unrepresentable width") + { + config.sectionWidth = std::numeric_limits::max(); + expectedCode = -7510; + } + SECTION("Zero columns") + { + config.sectionsPerRow = 0; + expectedCode = -7511; + } + SECTION("Too many columns") + { + config.sectionsPerRow = std::numeric_limits::max(); + expectedCode = -7511; + } + SECTION("Unrepresentable section count") + { + config.sectionCount = std::numeric_limits::max(); + expectedCode = -7510; + } + SECTION("Preparation error propagation") + { + config.grid.values = nullptr; + expectedCode = -7500; + } + SECTION("Nonfinite automatic MUD maximum") + { + values->initializeWithValue(std::numeric_limits::max()); + config.grid.units = ODFValueUnits::CountDensity; + expectedCode = -7512; + } + SECTION("Invalid manual endpoints") + { + config.scaleMode = ODFScaleMode::Manual; + config.manualMinimumMUD = 1; + config.manualMaximumMUD = 1; + expectedCode = -7512; + } + SECTION("Negative manual minimum") + { + config.scaleMode = ODFScaleMode::Manual; + config.manualMinimumMUD = -1; + expectedCode = -7512; + } + SECTION("Nonfinite manual maximum") + { + config.scaleMode = ODFScaleMode::Manual; + config.manualMaximumMUD = std::numeric_limits::infinity(); + expectedCode = -7512; + } + SECTION("Zero automatic maximum") + { + values->initializeWithValue(0.0); + expectedCode = -7512; + } + SECTION("One color control") + { + config.colorControlPoints = {0, 0, 0, 1}; + expectedCode = -7513; + } + SECTION("Incomplete control") + { + config.colorControlPoints.push_back(0); + expectedCode = -7513; + } + SECTION("Nonfinite control") + { + config.colorControlPoints[1] = std::numeric_limits::quiet_NaN(); + expectedCode = -7513; + } + SECTION("Repeated positions") + { + config.colorControlPoints[4] = 0; + expectedCode = -7513; + } + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == expectedCode); + REQUIRE_FALSE(result.errorMessage.empty()); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::RgbaPanelsAndChrome", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(32, "Planes", true); + values->initializeWithValue(0.0); + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < 2; phiIndex++) + { + const size_t flatIndex = (phi1Index * 2 + phiIndex) * 4 + 2; + values->setValue(flatIndex, 8.0); + } + } + const auto config = RenderConfiguration(values.get()); + const auto layout = ComputeODFSectionLayout(config); + REQUIRE(layout.panelHeight == 64); + REQUIRE(layout.panelSlotHeight == 118); + REQUIRE(layout.pageWidth == 534); + REQUIRE(layout.pageHeight == 270); + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE(result); + REQUIRE(result.image->getNumberOfComponents() == 4); + REQUIRE(result.image->getNumberOfTuples() == static_cast(result.width * result.height)); + REQUIRE(result.width == layout.pageWidth); + REQUIRE(result.height == layout.pageHeight); + REQUIRE(result.sectionAnglesDeg == std::vector{0, 90, 180, 270}); + REQUIRE(result.appliedMinimumMUD == Approx(0)); + REQUIRE(result.appliedMaximumMUD == Approx(4)); + RequirePixel(result, 64, 76, {0, 0, 255, 255}); + RequirePixel(result, 64, 194, {255, 0, 0, 255}); + RequirePixel(result, 480, 266, {255, 255, 255, 255}); + RequirePixel(result, 64, 108, {0, 0, 0, 255}); + RequirePixel(result, 362, 44, {255, 0, 0, 255}); + RequirePixel(result, 362, 107, {0, 0, 255, 255}); + REQUIRE(HasInk(result, 8, 8, 250, 12)); + REQUIRE(HasInk(result, 8, 26, 128, 12)); + REQUIRE(HasInk(result, 8, 121, 128, 12)); + REQUIRE(HasInk(result, 378, 54, 148, 84)); + for(size_t pixelIndex = 0; pixelIndex < result.image->getNumberOfTuples(); pixelIndex++) + { + REQUIRE(result.image->getValue(pixelIndex * 4 + 3) == 255); + } +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::SourceUnitsAndDenseChrome", "[EbsdLib][ODFSectionCompositor]") +{ + auto mudValues = DoubleArrayType::CreateArray(32, "MudValues", true); + mudValues->initializeWithValue(1.5); + auto densityValues = DoubleArrayType::CreateArray(32, "DensityValues", true); + const double stepRad = std::acos(-1.0) / 2.0; + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < 2; phiIndex++) + { + for(size_t phi2Index = 0; phi2Index < 4; phi2Index++) + { + const size_t flatIndex = (phi1Index * 2 + phiIndex) * 4 + phi2Index; + densityValues->setValue(flatIndex, 1.5 * stepRad * stepRad * stepRad * std::sin((phiIndex + 0.5) * stepRad) / (8.0 * std::acos(-1.0) * std::acos(-1.0))); + } + } + } + const auto mudConfig = RenderConfiguration(mudValues.get()); + auto densityConfig = RenderConfiguration(densityValues.get()); + densityConfig.grid.units = ODFValueUnits::CountDensity; + const auto layout = ComputeODFSectionLayout(mudConfig); + const auto mudResult = ODFSectionCompositor{}.generateCompositeImage(mudConfig); + const auto densityResult = ODFSectionCompositor{}.generateCompositeImage(densityConfig); + REQUIRE(mudResult); + REQUIRE(densityResult); + REQUIRE(mudResult.width == densityResult.width); + REQUIRE(mudResult.height == densityResult.height); + REQUIRE(mudResult.width == layout.pageWidth); + REQUIRE(mudResult.height == layout.pageHeight); + REQUIRE(mudResult.appliedMinimumMUD == Approx(densityResult.appliedMinimumMUD)); + REQUIRE(mudResult.appliedMaximumMUD == Approx(densityResult.appliedMaximumMUD)); + for(size_t sectionIndex = 0; sectionIndex < mudConfig.sectionCount; sectionIndex++) + { + const auto origin = GetODFSectionPanelOrigin(layout, sectionIndex); + for(int32_t row = 0; row < layout.panelHeight; row++) + { + for(int32_t column = 0; column < layout.panelWidth; column++) + { + const size_t offset = (static_cast(static_cast(origin[1]) + row) * mudResult.width + static_cast(origin[0]) + column) * 4; + for(size_t componentIndex = 0; componentIndex < 4; componentIndex++) + { + REQUIRE(mudResult.image->getValue(offset + componentIndex) == densityResult.image->getValue(offset + componentIndex)); + } + } + } + } + REQUIRE_FALSE(HasInk(mudResult, 400, 28, 80, 16)); + REQUIRE(HasInk(densityResult, 400, 28, 80, 16)); + REQUIRE(HasPixelDifference(mudResult, densityResult, 378, 84, 148, 17)); + REQUIRE(HasInk(mudResult, 65, 108, 7, 6)); + REQUIRE(HasInk(mudResult, 86, 125, 20, 13)); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::AppliedScaleAndColorInterpolation", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(32, "Uniform", true); + values->initializeWithValue(2.0); + auto config = RenderConfiguration(values.get()); + std::array expected = {128, 0, 128, 255}; + double expectedMinimum = 0; + double expectedMaximum = 4; + SECTION("Linear midpoint") + { + } + SECTION("Automatic uses displayed maximum") + { + config.scaleMode = ODFScaleMode::Automatic; + config.manualMinimumMUD = std::numeric_limits::quiet_NaN(); + expectedMaximum = 2; + expected = {255, 0, 0, 255}; + } + SECTION("Manual lower clipping") + { + config.manualMinimumMUD = 3; + expectedMinimum = 3; + expected = {0, 0, 255, 255}; + } + SECTION("Explicit zero grid remains valid with manual scale") + { + values->initializeWithValue(0.0); + expected = {0, 0, 255, 255}; + } + SECTION("Manual upper clipping") + { + config.manualMaximumMUD = 1; + expectedMaximum = 1; + expected = {255, 0, 0, 255}; + } + SECTION("Unequal control intervals") + { + config.colorControlPoints = {0, 0, 0, 1, 0.25f, 0, 1, 0, 1, 1, 0, 0}; + expected = {85, 170, 0, 255}; + } + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE(result); + REQUIRE(result.appliedMinimumMUD == Approx(expectedMinimum)); + REQUIRE(result.appliedMaximumMUD == Approx(expectedMaximum)); + RequirePixel(result, 64, 76, expected); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::AutomaticUsesCropAndInterpolation", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(32, "Cropped", true); + values->initializeWithValue(100.0); + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phi2Index = 0; phi2Index < 4; phi2Index++) + { + const size_t flatIndex = phi1Index * 8 + phi2Index; + values->setValue(flatIndex, static_cast(phi2Index + 1)); + } + } + auto config = RenderConfiguration(values.get()); + config.grid.laueOpsIndex = CrystalStructure::Hexagonal_High; + config.scaleMode = ODFScaleMode::Automatic; + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE(result); + REQUIRE(result.sectionAnglesDeg == std::vector{0, 15, 30, 45}); + REQUIRE(result.appliedMaximumMUD == Approx(2.5)); +} + +TEST_CASE("ebsdlib::ODFSectionCompositor::SpatialOrderAndUnusedSlot", "[EbsdLib][ODFSectionCompositor]") +{ + auto values = DoubleArrayType::CreateArray(32, "SpatialRamp", true); + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < 2; phiIndex++) + { + for(size_t phi2Index = 0; phi2Index < 4; phi2Index++) + { + const size_t flatIndex = (phi1Index * 2 + phiIndex) * 4 + phi2Index; + values->setValue(flatIndex, static_cast(phiIndex + phi1Index)); + } + } + } + auto config = RenderConfiguration(values.get()); + config.sectionCount = 3; + const auto result = ODFSectionCompositor{}.generateCompositeImage(config); + REQUIRE(result); + REQUIRE(result.sectionAnglesDeg == std::vector{0, 120, 240}); + RequirePixel(result, 49, 60, {0, 0, 255, 255}); + RequirePixel(result, 81, 60, {64, 0, 191, 255}); + RequirePixel(result, 145, 60, {191, 0, 64, 255}); + RequirePixel(result, 49, 92, {64, 0, 191, 255}); + RequirePixel(result, 145, 92, {255, 0, 0, 255}); + RequirePixel(result, 266, 194, {255, 255, 255, 255}); + // Annotations must not cover any pixels inside the blue source cell. + for(int32_t y = 46; y < 70; y++) + { + for(int32_t x = 35; x < 63; x++) + { + RequirePixel(result, x, y, {0, 0, 255, 255}); + } + } +} diff --git a/Source/Test/ODFSectionUtilitiesTest.cpp b/Source/Test/ODFSectionUtilitiesTest.cpp new file mode 100644 index 00000000..40a578db --- /dev/null +++ b/Source/Test/ODFSectionUtilitiesTest.cpp @@ -0,0 +1,311 @@ +/* ============================================================================ + * Copyright (c) 2009-2025 BlueQuartz Software, LLC + * + * Redistribution and use in source and binary forms, with or without modification, + * are permitted provided that the following conditions are met: + * + * Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * + * Redistributions in binary form must reproduce the above copyright notice, this + * list of conditions and the following disclaimer in the documentation and/or + * other materials provided with the distribution. + * + * Neither the name of BlueQuartz Software, the US Air Force, nor the names of its + * contributors may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR + * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER + * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, + * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE + * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + * + * The code contained herein was partially funded by the following contracts: + * United States Air Force Prime Contract FA8650-07-D-5800 + * United States Air Force Prime Contract FA8650-10-D-5210 + * United States Prime Contract Navy N00173-07-C-2068 + * + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ + +#include + +#include "EbsdLib/Math/EbsdLibMath.h" +#include "EbsdLib/Utilities/ODFSectionUtilities.h" + +#include +#include +#include + +using namespace ebsdlib; + +TEST_CASE("ebsdlib::ODFSectionUtilities::LaueLimitsAndAngles", "[EbsdLib][ODFSectionUtilities]") +{ + const auto limits = GetODFEulerPlotLimits(CrystalStructure::Hexagonal_High); + REQUIRE(limits.phi1MaxDeg == Approx(360.0)); + REQUIRE(limits.phiMaxDeg == Approx(90.0)); + REQUIRE(limits.phi2MaxDeg == Approx(60.0)); + REQUIRE(GenerateODFSectionAngles(CrystalStructure::Hexagonal_High, 6) == std::vector{0.0, 10.0, 20.0, 30.0, 40.0, 50.0}); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::CountDensityToMUD", "[EbsdLib][ODFSectionUtilities]") +{ + auto values = DoubleArrayType::CreateArray(32, "CountDensity", true); + values->initializeWithValue(1.0 / 32.0); + ODFGridView grid{values.get(), {4, 2, 4}, {0.0, 0.0, 0.0}, {90.0, 90.0, 90.0}, ODFValueUnits::CountDensity, CrystalStructure::Triclinic}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE(result); + const double step = constants::k_PiD / 2.0; + const double expected = (1.0 / 32.0) * 8.0 * constants::k_PiD * constants::k_PiD / (step * step * step * std::sin(constants::k_PiD / 4.0)); + REQUIRE(result.sections.mudValues.front() == Approx(expected)); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::Phi2FastestSourceLinearization", "[EbsdLib][ODFSectionUtilities]") +{ + constexpr size_t k_Phi1Count = 4; + constexpr size_t k_PhiCount = 2; + constexpr size_t k_Phi2Count = 4; + auto values = DoubleArrayType::CreateArray(k_Phi1Count * k_PhiCount * k_Phi2Count, "Phi2FastestSentinel", true); + for(size_t phi1Index = 0; phi1Index < k_Phi1Count; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < k_PhiCount; phiIndex++) + { + for(size_t phi2Index = 0; phi2Index < k_Phi2Count; phi2Index++) + { + const size_t flatIndex = (phi1Index * k_PhiCount + phiIndex) * k_Phi2Count + phi2Index; + values->setValue(flatIndex, 100.0 * static_cast(phi1Index) + 10.0 * static_cast(phiIndex) + static_cast(phi2Index)); + } + } + } + + const ODFGridView grid{values.get(), {k_Phi1Count, k_PhiCount, k_Phi2Count}, {0.0, 0.0, 0.0}, {90.0, 90.0, 90.0}, ODFValueUnits::MUD, CrystalStructure::Triclinic}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE(result); + REQUIRE(result.sections.mudValues == std::vector{1.5, 101.5, 201.5, 301.5, 11.5, 111.5, 211.5, 311.5, 0.5, 100.5, 200.5, 300.5, 10.5, 110.5, 210.5, 310.5, + 1.5, 101.5, 201.5, 301.5, 11.5, 111.5, 211.5, 311.5, 2.5, 102.5, 202.5, 302.5, 12.5, 112.5, 212.5, 312.5}); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::AllLaueLimits", "[EbsdLib][ODFSectionUtilities]") +{ + struct LimitsCase + { + uint32_t laueOpsIndex; + ODFEulerPlotLimits expected; + }; + + const std::array cases = {{{CrystalStructure::Hexagonal_High, {360.0, 90.0, 60.0}}, + {CrystalStructure::Cubic_High, {360.0, 90.0, 90.0}}, + {CrystalStructure::Hexagonal_Low, {360.0, 180.0, 60.0}}, + {CrystalStructure::Cubic_Low, {360.0, 90.0, 180.0}}, + {CrystalStructure::Triclinic, {360.0, 180.0, 360.0}}, + {CrystalStructure::Monoclinic, {360.0, 90.0, 360.0}}, + {CrystalStructure::OrthoRhombic, {360.0, 90.0, 180.0}}, + {CrystalStructure::Tetragonal_Low, {360.0, 180.0, 90.0}}, + {CrystalStructure::Tetragonal_High, {360.0, 90.0, 90.0}}, + {CrystalStructure::Trigonal_Low, {360.0, 180.0, 120.0}}, + {CrystalStructure::Trigonal_High, {360.0, 90.0, 120.0}}}}; + + for(const auto& testCase : cases) + { + const auto limits = GetODFEulerPlotLimits(testCase.laueOpsIndex); + REQUIRE(limits.phi1MaxDeg == Approx(testCase.expected.phi1MaxDeg)); + REQUIRE(limits.phiMaxDeg == Approx(testCase.expected.phiMaxDeg)); + REQUIRE(limits.phi2MaxDeg == Approx(testCase.expected.phi2MaxDeg)); + } +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::PeriodicPhi2Interpolation", "[EbsdLib][ODFSectionUtilities]") +{ + auto values = DoubleArrayType::CreateArray(32, "FourPlanes", true); + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < 2; phiIndex++) + { + for(size_t phi2Index = 0; phi2Index < 4; phi2Index++) + { + const size_t flatIndex = (phi1Index * 2 + phiIndex) * 4 + phi2Index; + values->setValue(flatIndex, static_cast(phi2Index + 1)); + } + } + } + + const ODFGridView grid{values.get(), {4, 2, 4}, {0.0, 0.0, 0.0}, {90.0, 90.0, 90.0}, ODFValueUnits::MUD, CrystalStructure::Triclinic}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE(result); + REQUIRE(result.sections.sectionAnglesDeg == std::vector{0.0, 90.0, 180.0, 270.0}); + REQUIRE(result.sections.phi1Count == 4); + REQUIRE(result.sections.phiCount == 2); + REQUIRE(result.sections.mudValues == + std::vector{2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5}); + REQUIRE(result.sections.mudValues.front() == Approx(2.5)); + REQUIRE(result.sections.maximumDisplayedMUD == Approx(3.5)); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::HexagonalHighCrop", "[EbsdLib][ODFSectionUtilities]") +{ + auto values = DoubleArrayType::CreateArray(32, "HexagonalHighCrop", true); + for(size_t phi1Index = 0; phi1Index < 4; phi1Index++) + { + for(size_t phiIndex = 0; phiIndex < 2; phiIndex++) + { + for(size_t phi2Index = 0; phi2Index < 4; phi2Index++) + { + const size_t flatIndex = (phi1Index * 2 + phiIndex) * 4 + phi2Index; + values->setValue(flatIndex, static_cast(phi1Index + 1 + phiIndex * 100)); + } + } + } + + const ODFGridView grid{values.get(), {4, 2, 4}, {0.0, 0.0, 0.0}, {90.0, 90.0, 90.0}, ODFValueUnits::MUD, CrystalStructure::Hexagonal_High}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE(result); + REQUIRE(result.sections.sectionAnglesDeg == std::vector{0.0, 15.0, 30.0, 45.0}); + REQUIRE(result.sections.phi1Count == 4); + REQUIRE(result.sections.phiCount == 1); + REQUIRE(result.sections.mudValues == std::vector{1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0, 1.0, 2.0, 3.0, 4.0}); + REQUIRE(result.sections.maximumDisplayedMUD == Approx(4.0)); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::EmptyDisplayedCrop", "[EbsdLib][ODFSectionUtilities]") +{ + auto values = DoubleArrayType::CreateArray(4, "CoarseGrid", true); + values->initializeWithValue(1.0); + const ODFGridView grid{values.get(), {2, 1, 2}, {0.0, 0.0, 0.0}, {180.0, 180.0, 180.0}, ODFValueUnits::MUD, CrystalStructure::Hexagonal_High}; + const auto result = PrepareODFSections(grid, 6); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + REQUIRE(result.errorMessage.find("(2, 1, 2)") != std::string::npos); + REQUIRE(result.errorMessage.find("180") != std::string::npos); + REQUIRE(result.errorMessage.find("PHI=0") != std::string::npos); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::Float32FractionalSpacing", "[EbsdLib][ODFSectionUtilities]") +{ + const size_t phiCount = GENERATE(size_t{7}, size_t{53}); + const double spacing = static_cast(180.0 / static_cast(phiCount)); + auto values = DoubleArrayType::CreateArray(4 * phiCount * phiCount * phiCount, "FractionalSpacing", true); + values->initializeWithValue(1.0); + ODFGridView grid{values.get(), {2 * phiCount, phiCount, 2 * phiCount}, {0.0, 0.0, 0.0}, {spacing, spacing, spacing}, ODFValueUnits::MUD, CrystalStructure::Hexagonal_High}; + const auto result = PrepareODFSections(grid, 6); + INFO(result.errorMessage); + REQUIRE(result); + REQUIRE_FALSE(result.sections.mudValues.empty()); + for(const double value : result.sections.mudValues) + { + REQUIRE(value == Approx(1.0)); + } + grid.spacingDeg = {spacing * 1.000001, spacing * 1.000001, spacing * 1.000001}; + const auto invalid = PrepareODFSections(grid, 6); + REQUIRE_FALSE(invalid); + REQUIRE(invalid.errorCode == -7502); +} + +TEST_CASE("ebsdlib::ODFSectionUtilities::InvalidInput", "[EbsdLib][ODFSectionUtilities]") +{ + auto values = DoubleArrayType::CreateArray(32, "ValidValues", true); + values->initializeWithValue(1.0); + const ODFGridView validGrid{values.get(), {4, 2, 4}, {0.0, 0.0, 0.0}, {90.0, 90.0, 90.0}, ODFValueUnits::MUD, CrystalStructure::Triclinic}; + + SECTION("Null values") + { + auto grid = validGrid; + grid.values = nullptr; + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7500); + } + + SECTION("Unsupported crystal code") + { + auto grid = validGrid; + grid.laueOpsIndex = CrystalStructure::UnknownCrystalStructure; + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7501); + } + + SECTION("Unsupported value units") + { + auto grid = validGrid; + grid.units = static_cast(2); + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + REQUIRE(result.errorMessage.find("(2)") != std::string::npos); + REQUIRE(result.errorMessage.find("MUD (0) or CountDensity (1)") != std::string::npos); + } + + SECTION("Transposed dimensions") + { + auto grid = validGrid; + grid.dimensions = {2, 4, 4}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + } + + SECTION("Nonzero origin") + { + auto grid = validGrid; + grid.originDeg = {0.0, 0.0, 1.0}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + } + + SECTION("Nonuniform spacing") + { + auto grid = validGrid; + grid.spacingDeg = {90.0, 45.0, 90.0}; + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7502); + } + + SECTION("Value count mismatch") + { + auto shortValues = DoubleArrayType::CreateArray(31, "ShortValues", true); + shortValues->initializeWithValue(1.0); + auto grid = validGrid; + grid.values = shortValues.get(); + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7503); + } + + SECTION("Negative value") + { + auto negativeValues = DoubleArrayType::CreateArray(32, "NegativeValues", true); + negativeValues->initializeWithValue(1.0); + negativeValues->setValue(7, -0.25); + auto grid = validGrid; + grid.values = negativeValues.get(); + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7504); + } + + SECTION("Not-a-number value") + { + auto nanValues = DoubleArrayType::CreateArray(32, "NaNValues", true); + nanValues->initializeWithValue(1.0); + nanValues->setValue(11, std::numeric_limits::quiet_NaN()); + auto grid = validGrid; + grid.values = nanValues.get(); + const auto result = PrepareODFSections(grid, 4); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7504); + } + + SECTION("One section") + { + const auto result = PrepareODFSections(validGrid, 1); + REQUIRE_FALSE(result); + REQUIRE(result.errorCode == -7505); + } +} diff --git a/Source/Test/PoleFigureCompositorTest.cpp b/Source/Test/PoleFigureCompositorTest.cpp index 89f06463..e42bbc43 100644 --- a/Source/Test/PoleFigureCompositorTest.cpp +++ b/Source/Test/PoleFigureCompositorTest.cpp @@ -279,7 +279,7 @@ TEST_CASE("ebsdlib::PoleFigureCompositorTest::All_Laue_Classes", "[EbsdLib][Pole { LaueOps::Pointer op = ops[opsIdx]; const std::string rpg = op->getRotationPointGroup(); - // Skip Triclinic (no FZ boundary) and duplicates (OrthoRhombicOps appears twice) + // Skip Triclinic (no FZ boundary) and rotation point groups already tested. if(rpg == "1" || tested.count(rpg) > 0) { continue; diff --git a/Source/Test/TextureTest.cpp b/Source/Test/TextureTest.cpp index d94a709f..e84a2c8a 100644 --- a/Source/Test/TextureTest.cpp +++ b/Source/Test/TextureTest.cpp @@ -34,8 +34,17 @@ * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ #include +#include +#include +#include +#include +#include +#include #include +#include +#include #include +#include #include #include "EbsdLib/Core/EbsdMacros.h" @@ -215,3 +224,355 @@ TEST_CASE("ebsdlib::TextureTest::DirectStructureMatrix", "[EbsdLib][DirectStruct } #endif + +namespace +{ +constexpr float k_Sigma3Angle = 60.0f * ebsdlib::constants::k_PiOver180F; +constexpr std::array k_Sigma3Axis = {1.0f, 1.0f, 1.0f}; + +std::vector createUniformCubicOdf() +{ + const Texture::ODFTableEntries entries; + return Texture::CalculateODFData>(entries, true); +} + +int sigma3MdfBin() +{ + CubicOps ops; + RodriguesDType rod = AxisAngleDType(k_Sigma3Axis[0], k_Sigma3Axis[1], k_Sigma3Axis[2], k_Sigma3Angle).toRodrigues(); + rod = ops.getMDFFZRod(rod); + return ops.getMisoBin(rod); +} + +std::vector calculateCubicMdf(const std::vector& inputWeights) +{ + std::vector angles(inputWeights.size(), k_Sigma3Angle); + std::vector axes(inputWeights.size() * 3); + for(size_t index = 0; index < inputWeights.size(); index++) + { + axes[index * 3] = k_Sigma3Axis[0]; + axes[index * 3 + 1] = k_Sigma3Axis[1]; + axes[index * 3 + 2] = k_Sigma3Axis[2]; + } + + std::vector weights = inputWeights; + const std::vector odf = createUniformCubicOdf(); + std::vector mdf; + Texture::CalculateMDFData(angles, axes, weights, odf, mdf, angles.size()); + return mdf; +} + +std::vector calculateSeededCubicMdf(const std::vector& inputWeights, uint64_t seed) +{ + std::vector angles(inputWeights.size(), k_Sigma3Angle); + std::vector axes(inputWeights.size() * 3); + for(size_t index = 0; index < inputWeights.size(); index++) + { + axes[index * 3] = k_Sigma3Axis[0]; + axes[index * 3 + 1] = k_Sigma3Axis[1]; + axes[index * 3 + 2] = k_Sigma3Axis[2]; + } + + std::vector weights = inputWeights; + const std::vector odf = createUniformCubicOdf(); + std::vector mdf; + std::mt19937_64 generator(seed); + Texture::CalculateMDFData(angles, axes, weights, odf, mdf, angles.size(), generator); + return mdf; +} + +void requireBitwiseEqual(const std::vector& first, const std::vector& second) +{ + REQUIRE(first.size() == second.size()); + for(size_t index = 0; index < first.size(); index++) + { + CAPTURE(index); + REQUIRE(std::bit_cast(first[index]) == std::bit_cast(second[index])); + } +} + +float sumMdf(const std::vector& mdf) +{ + return std::accumulate(mdf.cbegin(), mdf.cend(), 0.0f); +} +} // namespace + +TEST_CASE("ebsdlib::TextureTest::CalculateMDFData normalizes weighted targets", "[EbsdLib][TextureTest]") +{ + const int targetBin = sigma3MdfBin(); + + SECTION("One row reserves half of a cubic MDF") + { + const std::vector mdf = calculateCubicMdf({2916.0f}); + REQUIRE(mdf[targetBin] == Approx(0.5f).margin(1.0e-6f)); + REQUIRE(sumMdf(mdf) == Approx(1.0f).margin(1.0e-5f)); + } + + SECTION("An overflowing row clamps the reserved mass") + { + const std::vector mdf = calculateCubicMdf({500000.0f}); + REQUIRE(mdf[targetBin] == Approx(1.0f).margin(1.0e-6f)); + REQUIRE(sumMdf(mdf) == Approx(1.0f).margin(1.0e-5f)); + REQUIRE(std::none_of(mdf.cbegin(), mdf.cend(), [](float value) { return value < 0.0f; })); + } + + SECTION("Duplicate rows accumulate in the folded bin") + { + const std::vector mdf = calculateCubicMdf({1458.0f, 1458.0f}); + REQUIRE(mdf[targetBin] == Approx(0.5f).margin(1.0e-6f)); + REQUIRE(sumMdf(mdf) == Approx(1.0f).margin(1.0e-5f)); + } + + SECTION("Empty weights produce a random normalized MDF") + { + const std::vector mdf = calculateCubicMdf({}); + REQUIRE(sumMdf(mdf) == Approx(1.0f).margin(1.0e-5f)); + } +} + +TEST_CASE("ebsdlib::TextureTest::CalculateMDFData seeded generator is reproducible", "[EbsdLib][TextureTest]") +{ + constexpr uint64_t k_Seed = 0x5EED1234ULL; + + SECTION("Empty weights") + { + const std::vector first = calculateSeededCubicMdf({}, k_Seed); + const std::vector second = calculateSeededCubicMdf({}, k_Seed); + requireBitwiseEqual(first, second); + REQUIRE(sumMdf(first) == Approx(1.0f).margin(1.0e-5f)); + } + + SECTION("Weighted target") + { + const std::vector first = calculateSeededCubicMdf({2916.0f}, k_Seed); + const std::vector second = calculateSeededCubicMdf({2916.0f}, k_Seed); + requireBitwiseEqual(first, second); + REQUIRE(sumMdf(first) == Approx(1.0f).margin(1.0e-5f)); + } + + SECTION("Legacy overload") + { + const std::vector legacy = calculateCubicMdf({}); + REQUIRE_FALSE(legacy.empty()); + REQUIRE(sumMdf(legacy) == Approx(1.0f).margin(1.0e-5f)); + } + + SECTION("Seeded plot data") + { + std::vector firstMdf = calculateSeededCubicMdf({}, k_Seed); + std::vector secondMdf = firstMdf; + std::vector firstAngles; + std::vector firstFrequencies; + std::vector secondAngles; + std::vector secondFrequencies; + std::mt19937_64 firstGenerator(k_Seed); + std::mt19937_64 secondGenerator(k_Seed); + REQUIRE(StatsGen::GenMDFPlotData(firstMdf, firstAngles, firstFrequencies, 100000, firstGenerator) == 0); + REQUIRE(StatsGen::GenMDFPlotData(secondMdf, secondAngles, secondFrequencies, 100000, secondGenerator) == 0); + requireBitwiseEqual(firstAngles, secondAngles); + requireBitwiseEqual(firstFrequencies, secondFrequencies); + } +} + +TEST_CASE("ebsdlib::TextureTest::UniformOdfExcludesUnreachableBins", "[EbsdLib][TextureTest]") +{ + const auto odf = Texture::CalculateODFData>({}, true); + CHECK(odf.front() == 0.0); + CHECK(odf.back() == 0.0); + CHECK(std::accumulate(odf.begin(), odf.end(), 0.0) == Approx(1.0)); +} + +TEMPLATE_TEST_CASE("ebsdlib::TextureTest::UniformOdfSamplingHasNoClampFallbacks", "[EbsdLib][TextureTest][OdfInBall]", HexagonalOps, CubicOps, HexagonalLowOps, CubicLowOps, TriclinicOps, + MonoclinicOps, OrthoRhombicOps, TetragonalLowOps, TetragonalOps, TrigonalLowOps, TrigonalOps) +{ + TestType ops; + const auto odf = Texture::CalculateODFData>({}, true); + std::mt19937_64 generator(5489); + std::discrete_distribution selectBin(odf.begin(), odf.end()); + std::uniform_real_distribution offset(0.0, 1.0); + size_t nonFiniteCount = 0; + size_t above1799Count = 0; + size_t above180Count = 0; + ops.resetClampFallbackCount(); + constexpr size_t k_DrawCount = 200000; + for(size_t sampleIdx = 0; sampleIdx < k_DrawCount; ++sampleIdx) + { + const int bin = selectBin(generator); + double random[3] = {offset(generator), offset(generator), offset(generator)}; + const auto euler = ops.determineEulerAngles(random, bin); + if(!std::isfinite(euler[0]) || !std::isfinite(euler[1]) || !std::isfinite(euler[2])) + { + ++nonFiniteCount; + } + if constexpr(std::is_same_v) + { + const double angle = euler.toAxisAngle()[3] * 180.0 / constants::k_PiD; + nonFiniteCount += !std::isfinite(angle); + above1799Count += angle > 179.9; + above180Count += angle > 180.0; + } + } + std::cout << "ODF sampling " << ops.getNameOfClass() << ": draws=" << k_DrawCount << " clamp_fallbacks=" << ops.clampFallbackCount() << " non_finite=" << nonFiniteCount << std::endl; + if constexpr(std::is_same_v) + { + std::cout << "Triclinic angles: above_179.9=" << above1799Count << " above_180.0=" << above180Count << std::endl; + CHECK(above180Count == 0); + } + CHECK(nonFiniteCount == 0); + CHECK(ops.clampFallbackCount() == 0); +} + +namespace +{ +/** + * @struct OdfGridFixture + * @brief Independent grid dimensions and expected bin counts for one Laue class. + */ +struct OdfGridFixture +{ + const char* className; + std::array bins; + std::array halfWidths; + size_t nonzeroBins; + size_t partialBins; +}; + +// These independent grid fixtures retain the sampler's homochoric dimensions, +// including the 0.7f coefficient for the Monoclinic first axis. +const std::array k_OdfGridFixtures = {{ + {"HexagonalOps", {36, 36, 12}, {0.75366927563336727, 0.75366927563336727, 0.26060550051600867}, 15552, 0}, + {"CubicOps", {18, 18, 18}, {0.38867959478510306, 0.38867959478510306, 0.38867959478510306}, 5832, 0}, + {"HexagonalLowOps", {72, 72, 12}, {1.3306700394914688, 1.3306700394914688, 0.26060550051600867}, 49760, 3136}, + {"CubicLowOps", {36, 36, 36}, {0.75366927563336727, 0.75366927563336727, 0.75366927563336727}, 46656, 0}, + {"TriclinicOps", {72, 72, 72}, {1.3306700394914688, 1.3306700394914688, 1.3306700394914688}, 205704, 20336}, + {"MonoclinicOps", {72, 36, 72}, {1.3004169905036356, 0.75366927563336727, 1.3306700394914688}, 138808, 10136}, + {"OrthorhombicOps", {36, 36, 36}, {0.75366927563336727, 0.75366927563336727, 0.75366927563336727}, 46656, 0}, + {"TetragonalLowOps", {72, 72, 18}, {1.3306700394914688, 1.3306700394914688, 0.75366927563336727}, 68528, 6240}, + {"TetragonalOps", {36, 36, 18}, {0.75366927563336727, 0.75366927563336727, 0.38867959478510306}, 23328, 0}, + {"TrigonalLowOps", {72, 72, 24}, {1.3306700394914688, 1.3306700394914688, 0.26060550051600867}, 99384, 5976}, + {"TrigonalOps", {36, 36, 24}, {0.75366927563336727, 0.75366927563336727, 0.51410390002343753}, 31104, 0}, +}}; +} // namespace + +TEMPLATE_TEST_CASE("ebsdlib::TextureTest::OdfInBallFractionsMatchVolume", "[EbsdLib][TextureTest][OdfInBall]", HexagonalOps, CubicOps, HexagonalLowOps, CubicLowOps, TriclinicOps, MonoclinicOps, + OrthoRhombicOps, TetragonalLowOps, TetragonalOps, TrigonalLowOps, TrigonalOps) +{ + TestType ops; + const auto fixtureIter = std::find_if(k_OdfGridFixtures.begin(), k_OdfGridFixtures.end(), [&ops](const auto& fixture) { return fixture.className == ops.getNameOfClass(); }); + REQUIRE(fixtureIter != k_OdfGridFixtures.end()); + const auto& fixture = *fixtureIter; + REQUIRE(ops.getOdfNumBins() == fixture.bins); + CHECK(ops.odfBinInBallFraction(-1) == 0.0); + CHECK(ops.odfBinInBallFraction(static_cast(ops.getODFSize())) == 0.0); + CHECK_FALSE(ops.isOdfBinReachable(-1)); + CHECK_FALSE(ops.isOdfBinReachable(static_cast(ops.getODFSize()))); + + const auto odf = Texture::CalculateODFData>({}, true); + const auto floatOdf = Texture::CalculateODFData>({}, true); + REQUIRE(odf.size() == ops.getODFSize()); + CHECK(std::accumulate(odf.begin(), odf.end(), 0.0) == Approx(1.0).margin(1.0e-10)); + CHECK(std::accumulate(floatOdf.begin(), floatOdf.end(), 0.0) == Approx(1.0).margin(1.0e-6)); + + std::vector sampleBins; + for(size_t z : {size_t(0), fixture.bins[2] / 2, fixture.bins[2] - 1}) + { + for(size_t y : {size_t(0), fixture.bins[1] / 2, fixture.bins[1] - 1}) + { + for(size_t x : {size_t(0), fixture.bins[0] / 2, fixture.bins[0] - 1}) + { + sampleBins.push_back(static_cast((z * fixture.bins[1] + y) * fixture.bins[0] + x)); + } + } + } + + size_t zeroCount = 0; + size_t partialCount = 0; + size_t nonzeroCount = 0; + size_t invalidCount = 0; + double totalFraction = 0.0; + std::array boundarySamples = {}; + std::vector fractions(odf.size()); + for(size_t binIndex = 0; binIndex < odf.size(); ++binIndex) + { + const int bin = static_cast(binIndex); + const double fraction = ops.odfBinInBallFraction(bin); + fractions[binIndex] = fraction; + totalFraction += fraction; + zeroCount += fraction == 0.0; + partialCount += fraction > 0.0 && fraction < 1.0; + nonzeroCount += odf[binIndex] > 0.0; + invalidCount += !std::isfinite(fraction) || fraction < 0.0 || fraction > 1.0; + invalidCount += ops.isOdfBinReachable(bin) != (fraction > 0.0); + invalidCount += (odf[binIndex] == 0.0) != (fraction == 0.0); + if(fraction > 0.0 && fraction < 1.0) + { + const size_t bucket = static_cast(fraction * 10.0); + if(boundarySamples[bucket]++ < 4) + { + sampleBins.push_back(bin); + } + } + } + CHECK(invalidCount == 0); + CHECK(nonzeroCount == odf.size() - zeroCount); + CHECK(nonzeroCount == fixture.nonzeroBins); + CHECK(partialCount == fixture.partialBins); + std::cout << "ODF bins " << ops.getNameOfClass() << ": total=" << odf.size() << " zero=" << zeroCount << " partial=" << partialCount << " nonzero=" << nonzeroCount << std::endl; + + bool weightsMatch = true; + for(size_t binIndex = 0; binIndex < odf.size(); ++binIndex) + { + weightsMatch = weightsMatch && std::abs(odf[binIndex] - fractions[binIndex] / totalFraction) < 1.0e-14; + } + CHECK(weightsMatch); + + if constexpr(std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v || + std::is_same_v || std::is_same_v) + { + // The previous empty-entry calculation adds 1 to every bin, then divides by the bin count. + const std::vector previousOdf(ops.getODFSize(), 1.0 / static_cast(ops.getODFSize())); + const std::vector previousFloatOdf(ops.getODFSize(), 1.0f / static_cast(ops.getODFSize())); + CHECK(std::memcmp(odf.data(), previousOdf.data(), odf.size() * sizeof(double)) == 0); + CHECK(std::memcmp(floatOdf.data(), previousFloatOdf.data(), floatOdf.size() * sizeof(float)) == 0); + CHECK(zeroCount == 0); + CHECK(partialCount == 0); + } + + std::mt19937_64 generator(20260923); + std::uniform_real_distribution unit(0.0, 1.0); + constexpr size_t k_MonteCarloPoints = 2000; + const double radius = std::cbrt(3.0 * constants::k_PiD / 4.0); + for(int bin : sampleBins) + { + const size_t index = static_cast(bin); + const std::array cell = {index % fixture.bins[0], (index / fixture.bins[0]) % fixture.bins[1], index / (fixture.bins[0] * fixture.bins[1])}; + size_t insideCount = 0; + for(size_t sampleIdx = 0; sampleIdx < k_MonteCarloPoints; ++sampleIdx) + { + double radiusSquared = 0.0; + for(size_t axis = 0; axis < 3; ++axis) + { + const double coordinate = fixture.halfWidths[axis] * (2.0 * (static_cast(cell[axis]) + unit(generator)) / static_cast(fixture.bins[axis]) - 1.0); + radiusSquared += coordinate * coordinate; + } + insideCount += radiusSquared <= radius * radius; + } + const double measured = static_cast(insideCount) / static_cast(k_MonteCarloPoints); + INFO("class=" << ops.getNameOfClass() << " bin=" << bin << " Monte Carlo=" << measured); + CHECK(fractions[index] == Approx(measured).margin(0.05)); + } +} + +TEST_CASE("ebsdlib::TextureTest::HomochoricClampFallbackIsObservable", "[EbsdLib][TextureTest][OdfInBall]") +{ + TriclinicOps ops; + ops.resetClampFallbackCount(); + double random[3] = {0.5, 0.5, 0.5}; + const auto euler = ops.determineEulerAngles(random, 0); + CHECK(std::isfinite(euler[0])); + CHECK(std::isfinite(euler[1])); + CHECK(std::isfinite(euler[2])); + CHECK(ops.clampFallbackCount() == 1); + ops.resetClampFallbackCount(); + CHECK(ops.clampFallbackCount() == 0); +}