diff --git a/.github/workflows/clang-format.yml b/.github/workflows/format_pr.yml similarity index 70% rename from .github/workflows/clang-format.yml rename to .github/workflows/format_pr.yml index c7edb805..b59992a0 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/format_pr.yml @@ -8,7 +8,7 @@ on: jobs: clang_format_pr: - runs-on: ubuntu-22.04 + runs-on: ubuntu-24.04 steps: - uses: actions/checkout@v2 @@ -16,15 +16,20 @@ jobs: fetch-depth: 2 - name: Add Problem Matcher uses: ammaraskar/gcc-problem-matcher@a141586609e2a558729b99a8c574c048f7f56204 + - name: Install clang-format-19 + run: | + sudo apt-get update + sudo apt-get install clang-format-19 + clang-format-19 --version - name: Check Formatting id: check_format continue-on-error: true run: | - python3 scripts/clang_format.py --format-version 15 --commits HEAD^ HEAD + python3 scripts/clang_format.py --format-version 19 --commits HEAD^ HEAD - name: Apply Formatting if: steps.check_format.outcome != 'success' run: | - python3 scripts/clang_format.py --format-version 15 --modify --commits HEAD^ HEAD + python3 scripts/clang_format.py --format-version 19 --modify --commits HEAD^ HEAD - name: Add Suggestions if: steps.check_format.outcome != 'success' uses: reviewdog/action-suggester@v1 diff --git a/.github/workflows/format_push.yml b/.github/workflows/format_push.yml new file mode 100644 index 00000000..53502dc9 --- /dev/null +++ b/.github/workflows/format_push.yml @@ -0,0 +1,24 @@ +name: clang-format + +on: + push: + branches: + - develop + - master + +jobs: + clang_format: + runs-on: ubuntu-24.04 + + steps: + - uses: actions/checkout@v2 + - name: Add Problem Matcher + uses: ammaraskar/gcc-problem-matcher@a141586609e2a558729b99a8c574c048f7f56204 + - name: Install clang-format-19 + run: | + sudo apt-get update + sudo apt-get install clang-format-19 + clang-format-19 --version + - name: Check Formatting + run: | + python3 scripts/clang_format.py --format-version 19 diff --git a/Source/EbsdLib/Core/EbsdDataArray.cpp b/Source/EbsdLib/Core/EbsdDataArray.cpp index 01a26a0a..f66ad281 100644 --- a/Source/EbsdLib/Core/EbsdDataArray.cpp +++ b/Source/EbsdLib/Core/EbsdDataArray.cpp @@ -52,8 +52,8 @@ #define EBSD_BYTE_SWAP_32(x) _byteswap_ulong(x) #define EBSD_BYTE_SWAP_64(x) _byteswap_uint64(x) -#elif(defined(__clang__) && __has_builtin(__builtin_bswap32) && __has_builtin(__builtin_bswap64)) || (defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 3))) -#if(defined(__clang__) && __has_builtin(__builtin_bswap16)) || (defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 8))) +#elif (defined(__clang__) && __has_builtin(__builtin_bswap32) && __has_builtin(__builtin_bswap64)) || (defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 3))) +#if (defined(__clang__) && __has_builtin(__builtin_bswap16)) || (defined(__GNUC__) && (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 8))) #define EBSD_BYTE_SWAP_16(x) __builtin_bswap16(x) #else #define EBSD_BYTE_SWAP_16(x) __builtin_bswap32((x) << 16) diff --git a/Source/EbsdLib/Core/EbsdLibDLLExport.h b/Source/EbsdLib/Core/EbsdLibDLLExport.h index c1fe4bb8..df18dfb3 100644 --- a/Source/EbsdLib/Core/EbsdLibDLLExport.h +++ b/Source/EbsdLib/Core/EbsdLibDLLExport.h @@ -59,13 +59,13 @@ building on Windows. #if defined(EbsdLib_EXPORTS) /* Compiling the EbsdLib DLL/Dylib */ #if defined(_MSC_VER) /* MSVC Compiler Case */ #define EbsdLib_EXPORT __declspec(dllexport) -#elif(__GNUC__ >= 4) /* GCC 4.x has support for visibility options */ +#elif (__GNUC__ >= 4) /* GCC 4.x has support for visibility options */ #define EbsdLib_EXPORT __attribute__((visibility("default"))) #endif #else /* Importing the DLL into another project */ #if defined(_MSC_VER) /* MSVC Compiler Case */ #define EbsdLib_EXPORT __declspec(dllimport) -#elif(__GNUC__ >= 4) /* GCC 4.x has support for visibility options */ +#elif (__GNUC__ >= 4) /* GCC 4.x has support for visibility options */ #define EbsdLib_EXPORT __attribute__((visibility("default"))) #endif #endif diff --git a/Source/EbsdLib/IO/HKL/H5CtfImporter.cpp b/Source/EbsdLib/IO/HKL/H5CtfImporter.cpp index 5001c1c1..fb44f46c 100644 --- a/Source/EbsdLib/IO/HKL/H5CtfImporter.cpp +++ b/Source/EbsdLib/IO/HKL/H5CtfImporter.cpp @@ -332,7 +332,7 @@ int H5CtfImporter::writeSliceData(hid_t fileId, CtfReader& reader, int z, int ac if(nullptr == dataPtr) { assert(false); - } // We are going to crash here. I would rather crash than have bad data + } // We are going to crash here. I would rather crash than have bad data dataPtr = dataPtr + (actualSlice * dims[0]); // Put the pointer at the proper offset into the larger array WRITE_EBSD_DATA_ARRAY(reader, int, gid, name); } @@ -342,7 +342,7 @@ int H5CtfImporter::writeSliceData(hid_t fileId, CtfReader& reader, int z, int ac if(nullptr == dataPtr) { assert(false); - } // We are going to crash here. I would rather crash than have bad data + } // We are going to crash here. I would rather crash than have bad data dataPtr = dataPtr + (actualSlice * dims[0]); // Put the pointer at the proper offset into the larger array WRITE_EBSD_DATA_ARRAY(reader, float, gid, name); } diff --git a/Source/EbsdLib/LaueOps/CubicLowOps.cpp b/Source/EbsdLib/LaueOps/CubicLowOps.cpp index b29ddd9f..dbd80dfc 100644 --- a/Source/EbsdLib/LaueOps/CubicLowOps.cpp +++ b/Source/EbsdLib/LaueOps/CubicLowOps.cpp @@ -316,48 +316,26 @@ RodriguesDType CubicLowOps::getMDFFZRod(const RodriguesDType& inRod) const double n1 = ax[0]; double n2 = ax[1], n3 = ax[2], w = ax[3]; - double FZn1 = w, FZn2 = 0.0, FZn3 = 0.0, FZw = 0.0; + double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0; + double FZw = w; n1 = fabs(n1); n2 = fabs(n2); n3 = fabs(n3); - if(n1 > n2) + // The tetrahedral rotation group only provides the <111> 3-fold axes, so only + // cyclic permutations of the axis components are symmetry-equivalent. Rotate + // cyclically so the largest component is first; n2/n3 must NOT be sorted. + if(n2 >= n1 && n2 >= n3) { - if(n1 > n3) - { - FZn1 = n1; - if(n2 > n3) - { - FZn2 = n2, FZn3 = n3; - } - else - { - FZn2 = n3, FZn3 = n2; - } - } - else - { - FZn1 = n3, FZn2 = n1, FZn3 = n2; - } + FZn1 = n2, FZn2 = n3, FZn3 = n1; + } + else if(n3 >= n1 && n3 >= n2) + { + FZn1 = n3, FZn2 = n1, FZn3 = n2; } else { - if(n2 > n3) - { - FZn1 = n2; - if(n1 > n3) - { - FZn2 = n1, FZn3 = n3; - } - else - { - FZn2 = n3, FZn3 = n1; - } - } - else - { - FZn1 = n3, FZn2 = n2, FZn3 = n1; - } + FZn1 = n1, FZn2 = n2, FZn3 = n3; } return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); diff --git a/Source/EbsdLib/LaueOps/CubicOps.cpp b/Source/EbsdLib/LaueOps/CubicOps.cpp index eb1f20bb..c4624f1d 100644 --- a/Source/EbsdLib/LaueOps/CubicOps.cpp +++ b/Source/EbsdLib/LaueOps/CubicOps.cpp @@ -846,25 +846,25 @@ void CubicOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double double mag = loadx * loadx + loady * loady + loadz * loadz; mag = std::sqrt(mag); - theta1 = (loadx + loady + loadz) / (mag * 1.732f); + theta1 = (loadx + loady + loadz) / (mag * ebsdlib::constants::k_Sqrt3D); theta1 = std::fabs(theta1); - theta2 = (loadx + loady - loadz) / (mag * 1.732f); + theta2 = (loadx + loady - loadz) / (mag * ebsdlib::constants::k_Sqrt3D); theta2 = std::fabs(theta2); - theta3 = (loadx - loady + loadz) / (mag * 1.732f); + theta3 = (loadx - loady + loadz) / (mag * ebsdlib::constants::k_Sqrt3D); theta3 = std::fabs(theta3); - theta4 = (-loadx + loady + loadz) / (mag * 1.732f); + theta4 = (-loadx + loady + loadz) / (mag * ebsdlib::constants::k_Sqrt3D); theta4 = std::fabs(theta4); - lambda1 = (loadx + loady) / (mag * 1.414f); + lambda1 = (loadx + loady) / (mag * ebsdlib::constants::k_Sqrt2D); lambda1 = std::fabs(lambda1); - lambda2 = (loadx + loadz) / (mag * 1.414f); + lambda2 = (loadx + loadz) / (mag * ebsdlib::constants::k_Sqrt2D); lambda2 = std::fabs(lambda2); - lambda3 = (loadx - loady) / (mag * 1.414f); + lambda3 = (loadx - loady) / (mag * ebsdlib::constants::k_Sqrt2D); lambda3 = std::fabs(lambda3); - lambda4 = (loadx - loadz) / (mag * 1.414f); + lambda4 = (loadx - loadz) / (mag * ebsdlib::constants::k_Sqrt2D); lambda4 = std::fabs(lambda4); - lambda5 = (loady + loadz) / (mag * 1.414f); + lambda5 = (loady + loadz) / (mag * ebsdlib::constants::k_Sqrt2D); lambda5 = std::fabs(lambda5); - lambda6 = (loady - loadz) / (mag * 1.414f); + lambda6 = (loady - loadz) / (mag * ebsdlib::constants::k_Sqrt2D); lambda6 = std::fabs(lambda6); schmid1 = theta1 * lambda6; schmid2 = theta1 * lambda4; diff --git a/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp b/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp index bc24c4d8..90f00690 100644 --- a/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/HexagonalLowOps.cpp @@ -386,24 +386,14 @@ RodriguesDType HexagonalLowOps::getMDFFZRod(const RodriguesDType& inRod) const FZn1 = n1; FZn2 = n2; FZn3 = n3; - if(angle > 30.0) + // The 6/m rotation group is only the 6-fold about c (no in-plane 2-folds), so the + // axis azimuth folds into a plain 60 degree wedge with no mirror alternation { n1n2mag = std::sqrt(n1 * n1 + n2 * n2); - if(int(angle / 30) % 2 == 0) - { - FZw = angle - (30.0 * int(angle / 30.0)); - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); - } - else - { - FZw = angle - (30.0 * int(angle / 30.0)); - FZw = 30.0f - FZw; - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); - } + FZw = std::fmod(static_cast(angle), 60.0); + FZw = FZw * ebsdlib::constants::k_PiOver180D; + FZn1 = n1n2mag * std::cos(FZw); + FZn2 = n1n2mag * std::sin(FZw); } return AxisAngleDType(FZn1, FZn2, FZn3, w).toRodrigues(); @@ -532,6 +522,15 @@ int HexagonalLowOps::getOdfBin(const RodriguesDType& rod) const void HexagonalLowOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // Every output must be defined before the schmid comparison chain below, which only assigns to + // them when a candidate beats the incumbent. Without these, schmidfactor was READ uninitialized + // by the first `if(schmid1 > schmidfactor)`, and slipsys/angleComps were left untouched whenever + // no candidate won. + schmidfactor = 0.0; + slipsys = 0; + angleComps[0] = 0.0; + angleComps[1] = 0.0; + double theta1, theta2, theta3, theta4, theta5, theta6, theta7, theta8, theta9; double lambda1, lambda2, lambda3, lambda4, lambda5, lambda6, lambda7, lambda8, lambda9, lambda10; double schmid1, schmid2, schmid3, schmid4, schmid5, schmid6; diff --git a/Source/EbsdLib/LaueOps/HexagonalOps.cpp b/Source/EbsdLib/LaueOps/HexagonalOps.cpp index 94557f49..a9f7458e 100644 --- a/Source/EbsdLib/LaueOps/HexagonalOps.cpp +++ b/Source/EbsdLib/LaueOps/HexagonalOps.cpp @@ -86,8 +86,8 @@ namespace HexagonalHigh { constexpr std::array k_OdfNumBins = {36, 36, 12}; // Represents a 5Deg bin in homochoric space -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)), +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_PiD / 6.0) - std::sin(ebsdlib::constants::k_PiD / 6.0))), (1.0 / 3.0))}; 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)}; @@ -610,7 +610,14 @@ int HexagonalOps::getOdfBin(const RodriguesDType& rod) const void HexagonalOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // schmidfactor was already seeded here, but slipsys and angleComps were not: the comparison chain + // below only assigns to them when a candidate beats the incumbent, so a load direction for which + // every candidate is 0 left both outputs holding whatever the caller passed in. schmidfactor = 0.0; + slipsys = 0; + angleComps[0] = 0.0; + angleComps[1] = 0.0; + double theta1, theta2, theta3, theta4, theta5, theta6, theta7, theta8, theta9; double lambda1, lambda2, lambda3, lambda4, lambda5, lambda6, lambda7, lambda8, lambda9, lambda10; double schmid1, schmid2, schmid3, schmid4, schmid5, schmid6; diff --git a/Source/EbsdLib/LaueOps/MonoclinicOps.cpp b/Source/EbsdLib/LaueOps/MonoclinicOps.cpp index dd19774b..47746298 100644 --- a/Source/EbsdLib/LaueOps/MonoclinicOps.cpp +++ b/Source/EbsdLib/LaueOps/MonoclinicOps.cpp @@ -247,19 +247,31 @@ RodriguesDType MonoclinicOps::getODFFZRod(const RodriguesDType& rod) const // ----------------------------------------------------------------------------- RodriguesDType MonoclinicOps::getMDFFZRod(const RodriguesDType& inRod) const { - throw ebsdlib::method_not_implemented("MonoclinicOps::getMDFFZRod not implemented"); - // /// FIXME: Are we missing code for MonoclinicOps MDF FZ Rodrigues calculation? + RodriguesDType rod = LaueOps::_calcRodNearestOrigin(inRod); + AxisAngleDType ax = rod.toAxisAngle(); - // double w = 0.0, n1 = 0.0, n2 = 0.0, n3 = 0.0; - // double FZw = 0.0, FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0; - // - // OrientationType rod = LaueOps::_calcRodNearestOrigin(inRod); - // AxisAngleDType ax = rod.toAxisAngle(); - // n1 = ax[0]; - // n2 = ax[1], n3 = ax[2], w = ax[3]; - // - // - // return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); + // The 2/m rotation group is the single 2-fold about b (Y). Conjugation maps the + // misorientation axis (n1, n2, n3) to (-n1, n2, -n3) and switching symmetry + // negates it, so the fundamental sector is n2 >= 0 and n3 >= 0 with n1 free. + if(ax[1] < 0.0) + { + // combined 2-fold conjugation + switching: flips n2 only + ax[1] = -ax[1]; + } + if(ax[2] < 0.0) + { + // 2-fold conjugation: flips n1 and n3, leaves n2 + ax[0] = -ax[0]; + ax[2] = -ax[2]; + } + else if(ax[2] == 0.0 && ax[0] < 0.0) + { + // Equator tie-break: the 2-fold conjugation still relates (n1, n2, 0) and + // (-n1, n2, 0); pick the half-plane with n1 >= 0 + ax[0] = -ax[0]; + } + + return AxisAngleDType(ax[0], ax[1], ax[2], ax[3]).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -386,8 +398,13 @@ int MonoclinicOps::getOdfBin(const RodriguesDType& rod) const void MonoclinicOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void MonoclinicOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp b/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp index 14a95670..b8a75122 100644 --- a/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp +++ b/Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp @@ -258,20 +258,13 @@ RodriguesDType OrthoRhombicOps::getODFFZRod(const RodriguesDType& rod) const // ----------------------------------------------------------------------------- RodriguesDType OrthoRhombicOps::getMDFFZRod(const RodriguesDType& inRod) const { - throw ebsdlib::method_not_implemented("OrthoRhombicOps::getMDFFZRod not implemented"); - - double FZn1 = 0.0f, FZn2 = 0.0f, FZn3 = 0.0f, FZw = 0.0f; - RodriguesDType rod = _calcRodNearestOrigin(inRod); AxisAngleDType ax = rod.toAxisAngle(); - // double n1 = ax[0]; - // double n2 = ax[1]; - // double n3 = ax[2]; - // double w = ax[3]; - - /// FIXME: Are we missing code for OrthoRhombic MDF FZ Rodrigues calculation? - return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); + // The 222 rotation group's three orthogonal 2-folds combined with switching + // symmetry generate every sign combination of the misorientation axis, so the + // fundamental sector is the first octant. + return AxisAngleDType(std::fabs(ax[0]), std::fabs(ax[1]), std::fabs(ax[2]), ax[3]).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -397,8 +390,13 @@ int OrthoRhombicOps::getOdfBin(const RodriguesDType& rod) const void OrthoRhombicOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void OrthoRhombicOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp b/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp index f135f32c..5665aaa9 100644 --- a/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/TetragonalLowOps.cpp @@ -83,9 +83,14 @@ namespace TetragonalLow { constexpr std::array k_OdfNumBins = {72, 72, 18}; // Represents a 5Deg bin in homochoric space +// NOTE: the third dimension previously computed pow(0.75*(pi/4 - sin(pi/2)), 1/3), +// which is the cube root of a NEGATIVE number (NaN) -- a mismatched-angle bug +// inherited from legacy DREAM3D 6.5 OrientationLib. The homochoric half-width +// formula is 0.75*(theta - sin(theta)) with a single theta; for the 4-fold +// c-axis of 4/m that theta is pi/2, matching TetragonalOps (4/mmm). static const std::array k_OdfDimInitValue = {std::pow((0.75 * ((ebsdlib::constants::k_PiD)-std::sin((ebsdlib::constants::k_PiD)))), (1.0 / 3.0)), std::pow((0.75 * ((ebsdlib::constants::k_PiD)-std::sin((ebsdlib::constants::k_PiD)))), (1.0 / 3.0)), - std::pow((0.75 * ((ebsdlib::constants::k_PiOver4D)-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), k_OdfDimInitValue[1] / static_cast(k_OdfNumBins[1] / 2), k_OdfDimInitValue[2] / static_cast(k_OdfNumBins[2] / 2)}; @@ -259,17 +264,34 @@ RodriguesDType TetragonalLowOps::getODFFZRod(const RodriguesDType& rod) const // ----------------------------------------------------------------------------- RodriguesDType TetragonalLowOps::getMDFFZRod(const RodriguesDType& inRod) const { - double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0, FZw = 0.0; - RodriguesDType rod = _calcRodNearestOrigin(inRod); AxisAngleDType ax = rod.toAxisAngle(); - FZn1 = std::fabs(ax[0]); - FZn2 = std::fabs(ax[1]); - FZn3 = std::fabs(ax[2]); - FZw = ax[3]; + double n1 = ax[0]; + double n2 = ax[1]; + double n3 = ax[2]; + double w = ax[3]; + + // The 4/m rotation group is only the 4-fold about c (no in-plane 2-folds), so the + // axis folds to n3 >= 0 (via -C2z, which leaves the azimuth unchanged) and the + // azimuth into a plain 90 degree wedge with no mirror alternation + if(n3 < 0) + { + n3 = -n3; + } + double angle = 180.0 * std::atan2(n2, n1) * ebsdlib::constants::k_1OverPiD; + if(angle < 0) + { + angle = angle + 360.0; + } + { + const double n1n2mag = std::sqrt(n1 * n1 + n2 * n2); + const double azimuth = std::fmod(angle, 90.0) * ebsdlib::constants::k_PiOver180D; + n1 = n1n2mag * std::cos(azimuth); + n2 = n1n2mag * std::sin(azimuth); + } - return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); + return AxisAngleDType(n1, n2, n3, w).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -395,8 +417,13 @@ int TetragonalLowOps::getOdfBin(const RodriguesDType& rod) const void TetragonalLowOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void TetragonalLowOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/TetragonalOps.cpp b/Source/EbsdLib/LaueOps/TetragonalOps.cpp index f4296b0a..5a60b23e 100644 --- a/Source/EbsdLib/LaueOps/TetragonalOps.cpp +++ b/Source/EbsdLib/LaueOps/TetragonalOps.cpp @@ -297,6 +297,12 @@ RodriguesDType TetragonalOps::getMDFFZRod(const RodriguesDType& inRod) const FZn2 = std::fabs(ax[1]); FZn3 = std::fabs(ax[2]); FZw = ax[3]; + // The 422 rotation group's <110> 2-fold axes make (n1, n2) and (n2, n1) equivalent, + // so the octant folds further to the sector where n1 >= n2 + if(FZn2 > FZn1) + { + std::swap(FZn1, FZn2); + } return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); } @@ -424,8 +430,13 @@ int TetragonalOps::getOdfBin(const RodriguesDType& rod) const void TetragonalOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void TetragonalOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/TriclinicOps.cpp b/Source/EbsdLib/LaueOps/TriclinicOps.cpp index 46da46b8..c4b8abb6 100644 --- a/Source/EbsdLib/LaueOps/TriclinicOps.cpp +++ b/Source/EbsdLib/LaueOps/TriclinicOps.cpp @@ -244,14 +244,27 @@ RodriguesDType TriclinicOps::getODFFZRod(const RodriguesDType& rod) const // ----------------------------------------------------------------------------- RodriguesDType TriclinicOps::getMDFFZRod(const RodriguesDType& inRod) const { - throw ebsdlib::method_not_implemented("TriclinicOps::getMDFFZRod not implemented"); - RodriguesDType rod = LaueOps::_calcRodNearestOrigin(inRod); - AxisAngleDType ax = rod.toAxisAngle(); - /// FIXME: Are we missing code for TriclinicOps MDF FZ Rodrigues calculation? - return ax.toRodrigues(); + // The -1 Laue class has no rotational symmetry, so the only misorientation axis + // equivalence is switching symmetry (a misorientation and its inverse describe the + // same boundary and have negated axes). Fold to the upper hemisphere n3 >= 0. + if(ax[2] < 0.0) + { + ax[0] = -ax[0]; + ax[1] = -ax[1]; + ax[2] = -ax[2]; + } + else if(ax[2] == 0.0 && (ax[1] < 0.0 || (ax[1] == 0.0 && ax[0] < 0.0))) + { + // Equator tie-break: switching still relates (n1, n2, 0) and (-n1, -n2, 0); + // pick the half-plane with n2 > 0, or n1 >= 0 along the n2 == 0 line + ax[0] = -ax[0]; + ax[1] = -ax[1]; + } + + return AxisAngleDType(ax[0], ax[1], ax[2], ax[3]).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -377,8 +390,13 @@ int TriclinicOps::getOdfBin(const RodriguesDType& rod) const void TriclinicOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void TriclinicOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp b/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp index a3da9133..565bab74 100644 --- a/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp +++ b/Source/EbsdLib/LaueOps/TrigonalLowOps.cpp @@ -339,49 +339,44 @@ RodriguesDType TrigonalLowOps::getODFFZRod(const RodriguesDType& rod) const // ----------------------------------------------------------------------------- RodriguesDType TrigonalLowOps::getMDFFZRod(const RodriguesDType& inRod) const { - double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0, FZw = 0.0; - float n1n2mag = 0.0f; + double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0; + double n1n2mag = 0.0; RodriguesDType rod = _calcRodNearestOrigin(inRod); AxisAngleDType ax = rod.toAxisAngle(); - float denom = static_cast(std::sqrt(ax[0] * ax[0] + ax[1] * ax[1] + ax[2] * ax[2])); + double w = ax[3]; + + double denom = std::sqrt(ax[0] * ax[0] + ax[1] * ax[1] + ax[2] * ax[2]); ax[0] = ax[0] / denom; ax[1] = ax[1] / denom; - ax[1] = ax[2] / denom; + ax[2] = ax[2] / denom; if(ax[2] < 0) { ax[0] = -ax[0], ax[1] = -ax[1], ax[2] = -ax[2]; } - float angle = static_cast(180.0 * std::atan2(ax[1], ax[0]) * ebsdlib::constants::k_1OverPiD); + double angle = 180.0 * std::atan2(ax[1], ax[0]) * ebsdlib::constants::k_1OverPiD; if(angle < 0) { - angle = angle + 360.0f; + angle = angle + 360.0; } FZn1 = ax[0]; FZn2 = ax[1]; FZn3 = ax[2]; - if(angle > 60.0f) + // The -3 rotation group is only the 3-fold about c (no in-plane 2-folds), so the + // axis azimuth folds into a plain 120 degree wedge with no mirror alternation. + // On the equator (n3 == 0) switching symmetry acts within the plane and combines + // with the 3-fold into a 60 degree identification. { - n1n2mag = static_cast(std::sqrt(ax[0] * ax[0] + ax[1] * ax[1])); - if(int(angle / 60) % 2 == 0) - { - FZw = angle - (60.0f * int(angle / 60.0f)); - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); - } - else - { - FZw = angle - (60.0f * int(angle / 60.0f)); - FZw = 60.0f - FZw; - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); - } + const double sector = (ax[2] == 0.0) ? 60.0 : 120.0; + double azimuth = std::fmod(angle, sector); + n1n2mag = std::sqrt(ax[0] * ax[0] + ax[1] * ax[1]); + azimuth = azimuth * ebsdlib::constants::k_PiOver180D; + FZn1 = n1n2mag * std::cos(azimuth); + FZn2 = n1n2mag * std::sin(azimuth); } - return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); + return AxisAngleDType(FZn1, FZn2, FZn3, w).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -507,8 +502,13 @@ int TrigonalLowOps::getOdfBin(const RodriguesDType& rod) const void TrigonalLowOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void TrigonalLowOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/LaueOps/TrigonalOps.cpp b/Source/EbsdLib/LaueOps/TrigonalOps.cpp index 84d08be1..497f774c 100644 --- a/Source/EbsdLib/LaueOps/TrigonalOps.cpp +++ b/Source/EbsdLib/LaueOps/TrigonalOps.cpp @@ -360,8 +360,8 @@ RodriguesDType TrigonalOps::getODFFZRod(const RodriguesDType& rod) const RodriguesDType TrigonalOps::getMDFFZRod(const RodriguesDType& inRod) const { double w = 0.0, n1 = 0.0, n2 = 0.0, n3 = 0.0; - double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0, FZw = 0.0; - double n1n2mag = 0.0f; + double FZn1 = 0.0, FZn2 = 0.0, FZn3 = 0.0; + double n1n2mag = 0.0; RodriguesDType rod = _calcRodNearestOrigin(inRod); @@ -386,27 +386,33 @@ RodriguesDType TrigonalOps::getMDFFZRod(const RodriguesDType& inRod) const FZn1 = n1; FZn2 = n2; FZn3 = n3; - if(angle > 60.0f) + // The 32 rotation group's in-plane 2-folds lie at azimuths 0/60/120 degrees, which + // places the mirror lines of the misorientation axis group at 30 + 60k degrees. The + // fundamental azimuth sector is therefore [30, 90] degrees. { - n1n2mag = std::sqrt(n1 * n1 + n2 * n2); - if(int(angle / 60) % 2 == 0) + double azimuth = angle - (120.0 * static_cast(angle / 120.0)); + if(azimuth < 30.0) + { + azimuth = 60.0 - azimuth; + } + else if(azimuth > 90.0) { - FZw = angle - (60.0f * int(angle / 60.0f)); - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); + azimuth = 180.0 - azimuth; } - else + if(n3 == 0.0 && azimuth > 60.0) { - FZw = angle - (60.0f * int(angle / 60.0f)); - FZw = 60.0f - FZw; - FZw = FZw * ebsdlib::constants::k_PiOver180D; - FZn1 = n1n2mag * std::cos(FZw); - FZn2 = n1n2mag * std::sin(FZw); + // Equator tie-break: on the equator the 2-fold conjugations act within the + // plane, adding mirror lines every 30 degrees; the sector shrinks to [30, 60] + azimuth = 120.0 - azimuth; } + n1n2mag = std::sqrt(n1 * n1 + n2 * n2); + azimuth = azimuth * ebsdlib::constants::k_PiOver180D; + FZn1 = n1n2mag * std::cos(azimuth); + FZn2 = n1n2mag * std::sin(azimuth); } - return AxisAngleDType(FZn1, FZn2, FZn3, FZw).toRodrigues(); + // The misorientation angle w is unchanged by the axis fold and must be preserved + return AxisAngleDType(FZn1, FZn2, FZn3, w).toRodrigues(); } // ----------------------------------------------------------------------------- @@ -531,8 +537,13 @@ int TrigonalOps::getOdfBin(const RodriguesDType& rod) const void TrigonalOps::getSchmidFactorAndSS(double load[3], double& schmidfactor, double angleComps[2], int& slipsys) const { + // No slip systems are enumerated for this Laue class. Zero EVERY output, angleComps + // included: leaving them untouched handed the caller back whatever it passed in, which for + // a caller that reuses one angleComps buffer across a loop is the PREVIOUS entry's angles. schmidfactor = 0; slipsys = 0; + angleComps[0] = 0; + angleComps[1] = 0; } void TrigonalOps::getSchmidFactorAndSS(double load[3], double plane[3], double direction[3], double& schmidfactor, double angleComps[2], int& slipsys) const diff --git a/Source/EbsdLib/Texture/MisorientationKDE.cpp b/Source/EbsdLib/Texture/MisorientationKDE.cpp new file mode 100644 index 00000000..76210841 --- /dev/null +++ b/Source/EbsdLib/Texture/MisorientationKDE.cpp @@ -0,0 +1,242 @@ +/* ============================================================================ + * 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 "MisorientationKDE.h" + +#include "EbsdLib/Orientation/Homochoric.hpp" +#include "EbsdLib/Orientation/Rodrigues.hpp" +#include "EbsdLib/Texture/RandomAngleDistribution.h" + +#include +#include +#include + +namespace ebsdlib +{ +MisorientationKDE::MisorientationKDE(LaueOps::Pointer ops, uint32_t crystalStructure, double halfwidthRadians) +: m_Ops(std::move(ops)) +, m_CrystalStructure(crystalStructure) +, m_Kernel(halfwidthRadians) +, m_BinWeights(m_Ops->getMDFSize(), 0.0) +, m_BinQuatSum(m_Ops->getMDFSize(), std::array{0.0, 0.0, 0.0, 0.0}) +{ + size_t numSymOps = m_Ops->getNumSymOps(); + m_SymQuats.reserve(numSymOps); + for(size_t i = 0; i < numSymOps; i++) + { + m_SymQuats.push_back(m_Ops->getQuatSymOp(i)); + } +} + +void MisorientationKDE::addMisorientation(const QuatD& misoQuat, double weight) +{ + RodriguesDType rod = m_Ops->getMDFFZRod(misoQuat.toRodrigues()); + int binIndex = m_Ops->getMisoBin(rod); + const size_t bin = static_cast(binIndex); + m_BinWeights[bin] += weight; + m_TotalWeight += weight; + + // Accumulate the weighted circular mean of the FZ misorientation quaternion within its bin. + // All misorientations sharing a bin lie within one ~5-degree cell, so a sign-aligned linear + // sum (aligned to the bin's running accumulator, or to itself when the bin is first seen) and + // a final renormalization is an accurate mean on that small patch of SO(3). Using this mean as + // the KDE center (in finalize) instead of the geometric bin center removes the bin-quantization + // bias without storing the individual misorientations (memory stays bounded by getMDFSize). + const QuatD fzQuat = rod.toQuaternion(); + std::array& acc = m_BinQuatSum[bin]; + const double alignDot = fzQuat.x() * acc[0] + fzQuat.y() * acc[1] + fzQuat.z() * acc[2] + fzQuat.w() * acc[3]; + const double sign = (alignDot < 0.0) ? -1.0 : 1.0; + acc[0] += weight * sign * fzQuat.x(); + acc[1] += weight * sign * fzQuat.y(); + acc[2] += weight * sign * fzQuat.z(); + acc[3] += weight * sign * fzQuat.w(); +} + +void MisorientationKDE::finalize() +{ + m_Centers.clear(); + if(m_TotalWeight <= 0.0) + { + return; + } + for(size_t binIndex = 0; binIndex < m_BinWeights.size(); binIndex++) + { + if(m_BinWeights[binIndex] > 0.0) + { + // Use the bin's weighted circular-mean quaternion (accumulated in addMisorientation) as the + // representative center rather than the geometric bin center. This eliminates the ~5-degree + // MDF-bin quantization bias that otherwise shifts the extracted angle-distribution curve. + const std::array& acc = m_BinQuatSum[binIndex]; + const double norm = std::sqrt(acc[0] * acc[0] + acc[1] * acc[1] + acc[2] * acc[2] + acc[3] * acc[3]); + QuatD quat = (norm > 0.0) ? QuatD(acc[0] / norm, acc[1] / norm, acc[2] / norm, acc[3] / norm) : binCenter(static_cast(binIndex)); + m_Centers.push_back({quat, quat.conjugate(), m_BinWeights[binIndex] / m_TotalWeight}); + } + } +} + +double MisorientationKDE::totalWeight() const +{ + return m_TotalWeight; +} + +QuatD MisorientationKDE::binCenter(int binIndex) const +{ + double center[3] = {0.5, 0.5, 0.5}; + RodriguesDType rod = m_Ops->determineRodriguesVector(center, binIndex); + return rod.toQuaternion(); +} + +double MisorientationKDE::evaluate(const QuatD& query) const +{ + const double cutoffCos = std::cos(m_Kernel.cutoffAngle() / 2.0); + const size_t numSymOps = m_SymQuats.size(); + + // Symmetrize the query once: s_i * q * s_j for all crystal-symmetry pairs. + std::vector symQueries; + symQueries.reserve(numSymOps * numSymOps); + for(size_t i = 0; i < numSymOps; i++) + { + QuatD left = m_SymQuats[i] * query; + for(size_t j = 0; j < numSymOps; j++) + { + symQueries.push_back(left * m_SymQuats[j]); + } + } + + double density = 0.0; + for(const Center& center : m_Centers) + { + double kernelSum = 0.0; + for(const QuatD& symQuery : symQueries) + { + double dotForward = std::fabs(symQuery.dotProduct(center.Quat)); + if(dotForward >= cutoffCos) + { + kernelSum += m_Kernel.evaluate(dotForward); + } + double dotInverse = std::fabs(symQuery.dotProduct(center.QuatInverse)); + if(dotInverse >= cutoffCos) + { + kernelSum += m_Kernel.evaluate(dotInverse); + } + } + // Sum of the forward K(g, c) and grain-exchange K(g, inv(c)) kernels, averaged + // over the |CS| x |CS| crystal-symmetry pairs, then halved. The forward and inverse + // terms together span the grain-exchange-extended symmetry orbit (2 * |CS| * |CS| + // elements), so the correct mean-1 normalization divides by 2 * |CS|^2, i.e. the + // 0.5 antipodal factor. This was pinned by a direct numerical cross-check against + // MTEX 6.1.0 (calcDensity 'exact'): with exact (un-gridified) centers our density + // equals MTEX's mdf to a constant ratio of exactly 2.000 across all misorientation + // angles, and MTEX's mdf has mean(mdf) == 1 over SO(3). Dropping this factor gives a + // mean of ~2 (and a triclinic modal peak of K(0) instead of K(0)/2). See + // MisorientationKDETest.cpp::CubicAngleCurveVsMTEX for the cross-check. + density += center.Weight * kernelSum / static_cast(2 * numSymOps * numSymOps); + } + return density; +} + +std::vector MisorientationKDE::evaluateAtBinCenters() const +{ + size_t mdfSize = m_Ops->getMDFSize(); + std::vector densities(mdfSize, 0.0); + for(size_t binIndex = 0; binIndex < mdfSize; binIndex++) + { + densities[binIndex] = evaluate(binCenter(static_cast(binIndex))); + } + return densities; +} + +MisorientationKDE::AngleCurve MisorientationKDE::computeAngleCurve(size_t numPoints) const +{ + AngleCurve curve; + const uint32_t structure = m_CrystalStructure; + const double maxAngle = random_angle_distribution::MaxMisorientationAngle(structure); + + curve.Angles.resize(numPoints); + for(size_t i = 0; i < numPoints; i++) + { + curve.Angles[i] = maxAngle * static_cast(i) / static_cast(numPoints - 1); + } + curve.RandomDensity = random_angle_distribution::Compute(structure, curve.Angles); + curve.Density = curve.RandomDensity; // start from the uniform reference, MTEX-style + + const double resolution = 0.5 * constants::k_DegToRadD; // MTEX default 'resolution' + const double gridScale = 2.0 * static_cast(m_SymQuats.size()); // full-sphere grid vs MTEX sector grid + const size_t maxAxes = 20000; + constexpr double k_GoldenAngle = 2.399963229728653; + + for(size_t i = 0; i < numPoints; i++) + { + const double omega = curve.Angles[i]; + const double sinHalf = std::sin(omega / 2.0); + const double cosHalf = std::cos(omega / 2.0); + size_t numAxes = static_cast(std::lround(gridScale * (4.0 / 3.0) * sinHalf * sinHalf / (resolution * resolution))); + numAxes = std::clamp(numAxes, 1, maxAxes); + + double sum = 0.0; + size_t accepted = 0; + for(size_t a = 0; a < numAxes; a++) + { + // Fibonacci sphere point a of numAxes. + const double z = 1.0 - 2.0 * (static_cast(a) + 0.5) / static_cast(numAxes); + const double r = std::sqrt(std::max(0.0, 1.0 - z * z)); + const double phi = static_cast(a) * k_GoldenAngle; + const double axisX = r * std::cos(phi); + const double axisY = r * std::sin(phi); + const double axisZ = z; + + QuatD q(axisX * sinHalf, axisY * sinHalf, axisZ * sinHalf, cosHalf); + // Keep only MDF-fundamental-zone representatives: fold with getMDFFZRod and + // compare in homochoric space; each misorientation class is counted once. + RodriguesDType rod = q.toRodrigues(); + HomochoricDType hoOriginal = rod.toHomochoric(); + HomochoricDType hoFolded = m_Ops->getMDFFZRod(rod).toHomochoric(); + const double tolerance = 1.0e-6; + if(std::fabs(hoOriginal[0] - hoFolded[0]) > tolerance || std::fabs(hoOriginal[1] - hoFolded[1]) > tolerance || std::fabs(hoOriginal[2] - hoFolded[2]) > tolerance) + { + continue; + } + sum += evaluate(q); + accepted++; + } + if(accepted > 0) + { + curve.Density[i] *= std::max(0.0, sum / static_cast(accepted)); + } + // MTEX leaves density(i) at the uniform value when the slice grid is empty. + } + return curve; +} +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/MisorientationKDE.h b/Source/EbsdLib/Texture/MisorientationKDE.h new file mode 100644 index 00000000..244134dc --- /dev/null +++ b/Source/EbsdLib/Texture/MisorientationKDE.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/EbsdLib.h" +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/Orientation/Quaternion.hpp" +#include "EbsdLib/Texture/SO3DeLaValleePoussinKernel.h" + +#include +#include +#include + +namespace ebsdlib +{ +/** + * @brief Misorientation kernel density estimator on the MDF (misorientation) fundamental-zone bin grid. + * + * Accumulates weighted misorientations into the Laue-class MDF-FZ bins (via LaueOps::getMDFFZRod + + * getMisoBin), then evaluates a symmetrized De la Vallee Poussin kernel density at an arbitrary + * misorientation. The density is: + * + * f(g) = sum_bins w_bin * ( K(g, c_bin) + K(g, c_bin^-1) ) / ( 2 * |CS|^2 ) + * + * where each K is summed over the |CS| x |CS| crystal-symmetry pairs (s_i * g * s_j), and the + * grain-exchange (antipodal) inverse term enforces f(g) == f(g^-1). The forward and inverse terms + * together span the grain-exchange-extended symmetry orbit (2 * |CS|^2 elements), so the mean-1 + * normalization divides by 2 * |CS|^2. The kernel psi integrates to 1 over SO(3); with weights + * normalized to sum 1 the density is then a normalized MDF (mean == 1 over SO(3), matching MTEX's + * mean(mdf) == 1) whose triclinic (|CS| == 1) modal peak height is K(0) / 2. The absolute scale was + * pinned by a direct numerical cross-check against MTEX 6.1.0 (see computeAngleCurve()). + * + * Usage: construct, addMisorientation() for every observation, finalize() once, then evaluate(). + */ +class EbsdLib_EXPORT MisorientationKDE +{ +public: + /** + * @brief Misorientation-angle-distribution curve extracted from the KDE. + * Angles are in radians (0 .. MaxMisorientationAngle(structure)); Density is the MDF + * angle distribution; RandomDensity is the uniform (random) reference distribution. + */ + struct AngleCurve + { + std::vector Angles; + std::vector Density; + std::vector RandomDensity; + }; + + /** + * @brief Constructor. + * @param ops Laue-class symmetry operators for the MDF fundamental zone. + * @param crystalStructure EbsdLib crystal-structure index of ops. LaueOps has no reverse lookup; + * it is stored for Task 4's computeAngleCurve() (the Mackenzie reference) and is not + * consumed by this class. + * @param halfwidthRadians De la Vallee Poussin kernel halfwidth in radians. + */ + MisorientationKDE(LaueOps::Pointer ops, uint32_t crystalStructure, double halfwidthRadians); + + /** + * @brief Accumulate a weighted misorientation into its MDF-FZ bin. + * @param misoQuat Misorientation quaternion. + * @param weight Non-negative weight (need not be normalized). + */ + void addMisorientation(const QuatD& misoQuat, double weight); + + /** + * @brief Normalize accumulated weights to sum 1 and build the center list. Call once, after all adds. + */ + void finalize(); + + /** + * @brief Sum of all weights passed to addMisorientation() (before normalization). + */ + double totalWeight() const; + + /** + * @brief Density at an arbitrary misorientation quaternion; valid after finalize(). + */ + double evaluate(const QuatD& query) const; + + /** + * @brief MDF-FZ-folded bin-center quaternion for the given miso bin index. + */ + QuatD binCenter(int binIndex) const; + + /** + * @brief Serial convenience: evaluate() at every bin center; size getMDFSize(). + */ + std::vector evaluateAtBinCenters() const; + + /** + * @brief Misorientation-angle-distribution curve: a port of MTEX + * SO3Fun/@SO3Fun/calcAngleDistribution.m. For each of numPoints angles omega in + * [0, MaxMisorientationAngle(structure)] the density is the uniform-reference value + * (random_angle_distribution::Compute) multiplied by the mean of evaluate() over an + * axis grid on the omega-sphere. Two deliberate deviations from MTEX (validated by the + * MTEX numerical cross-check in the unit test, tolerance epsilon(0.20)+margin(0.10)): + * (1) axes are a Fibonacci full-sphere sampling filtered to MDF-FZ membership (via the + * audited getMDFFZRod folds, each misorientation class counted exactly once) + * instead of MTEX's fundamental-sector grid with a one-sided FZ check; + * (2) the axis-grid count is scaled by 2*|CS| so the post-filter (FZ-only) axis density + * matches MTEX's sector-grid density. + */ + AngleCurve computeAngleCurve(size_t numPoints) const; + +private: + struct Center + { + QuatD Quat; + QuatD QuatInverse; + double Weight; + }; + + LaueOps::Pointer m_Ops; + uint32_t m_CrystalStructure; + SO3DeLaValleePoussinKernel m_Kernel; + std::vector m_SymQuats; + std::vector m_BinWeights; + // Per-bin weight-weighted running sum of the (sign-aligned) fundamental-zone misorientation + // quaternions {x,y,z,w}. finalize() normalizes each non-empty bin's sum to obtain the bin's + // representative center. This is the weighted circular mean of the misorientations that fell + // in the bin, which is far closer to the true data than the geometric bin center and removes + // the ~5-degree MDF-bin quantization bias from the extracted angle-distribution curve. + std::vector> m_BinQuatSum; + double m_TotalWeight = 0.0; + std::vector
m_Centers; +}; +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/RandomAngleDistribution.cpp b/Source/EbsdLib/Texture/RandomAngleDistribution.cpp new file mode 100644 index 00000000..a83ac0c4 --- /dev/null +++ b/Source/EbsdLib/Texture/RandomAngleDistribution.cpp @@ -0,0 +1,293 @@ +/* ============================================================================ + * 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 "RandomAngleDistribution.h" + +#include "EbsdLib/Core/EbsdLibConstants.h" +#include "EbsdLib/Math/EbsdLibMath.h" + +#include +#include + +namespace +{ +// ----------------------------------------------------------------------------- +// Direct ports of the helper functions in MTEX 6.1.0's +// geometry/@symmetry/calcAngleDistribution.m (lines 231-255). + +// the area of the spherical triangle; alpha, beta, gamma are angles between vertices +double C(double alpha, double beta, double gamma) +{ + return std::acos((std::cos(gamma) - std::cos(alpha) * std::cos(beta)) / (std::sin(alpha) * std::sin(beta))); +} + +// the area of a spherical cap +double S1(double rho) +{ + return 2.0 * ebsdlib::constants::k_PiD * (1.0 - std::cos(rho)); +} + +// area of the intersection of two spherical caps; rho1, rho2 are radii of the +// caps, xi is the distance between the centers of the caps +double S2(double rho1, double rho2, double xi) +{ + return 2.0 * (ebsdlib::constants::k_PiD - C(rho1, rho2, xi) - std::cos(rho1) * C(xi, rho1, rho2) - std::cos(rho2) * C(rho2, xi, rho1)); +} + +// Cn branch (calcAngleDistribution.m lines 58-63): nfold is the crystal's +// rotational symmetry order about its principal axis. +double ChiCn(double rmag, uint32_t nfold) +{ + const double xhn = std::tan(ebsdlib::constants::k_PiD / 2.0 / static_cast(nfold)); + double xchi = 1.0; + if(rmag > xhn) + { + xchi = xhn / rmag; + } + return xchi; +} + +// Dnh branch (calcAngleDistribution.m lines 64-83), built on top of the Cn +// first region. +double ChiDnh(double rmag, uint32_t nfold) +{ + const double xhn = std::tan(ebsdlib::constants::k_PiD / 2.0 / static_cast(nfold)); + double xchi = ChiCn(rmag, nfold); + + if(rmag > 1.0) + { + xchi += static_cast(nfold) * (1.0 / rmag - 1.0); + } + + const double xedge = std::sqrt(1.0 + xhn * xhn); + if(rmag > xedge) + { + const double alpha1 = std::acos(xhn / rmag); + const double alpha2 = std::acos(1.0 / rmag); + const double xs21 = S2(alpha1, alpha2, ebsdlib::constants::k_PiD / 2.0); + const double xs22 = S2(alpha2, alpha2, ebsdlib::constants::k_PiD / static_cast(nfold)); + xchi += static_cast(nfold) * xs21 / ebsdlib::constants::k_PiD + static_cast(nfold) * xs22 / (2.0 * ebsdlib::constants::k_PiD); + } + + return xchi; +} + +// m-3 branch (calcAngleDistribution.m lines 85-96). +double ChiM3(double rmag) +{ + double xchi = 1.0; + + // first region + const double xh3 = std::sqrt(3.0) / 3.0; + if(rmag > xh3) + { + xchi = 4.0 * xh3 / rmag - 3.0; + } + + // second region + const double xedge = std::sqrt(2.0) / 2.0; + if(rmag > xedge) + { + const double alpha = std::acos(xh3 / rmag); + xchi += 3.0 * S2(alpha, alpha, std::acos(1.0 / 3.0)) / ebsdlib::constants::k_PiD; + } + + return xchi; +} + +// m-3m branch (calcAngleDistribution.m lines 98-117). +double ChiM3m(double rmag) +{ + double xchi = 1.0; + + // first region -> four fold axis active + const double xh4 = std::sqrt(2.0) - 1.0; + if(rmag > xh4) + { + xchi = 3.0 * xh4 / rmag - 2.0; + } + + // second region -> three fold axis active + const double xh3 = std::sqrt(3.0) / 3.0; + if(rmag > xh3) + { + xchi += 4.0 * (xh3 / rmag - 1.0); + } + + // third region + const double xedge = 2.0 - std::sqrt(2.0); + if(rmag > xedge) + { + const double alpha1 = std::acos(xh4 / rmag); + const double alpha2 = std::acos(xh3 / rmag); + const double s12 = S2(alpha1, alpha1, ebsdlib::constants::k_PiD / 2.0); + const double s24 = S2(alpha1, alpha2, std::acos(xh3)); + xchi += 3.0 * s12 / ebsdlib::constants::k_PiD + 6.0 * s24 / ebsdlib::constants::k_PiD; + } + + return xchi; +} +} // namespace + +namespace ebsdlib +{ +namespace random_angle_distribution +{ +double MaxMisorientationAngle(uint32_t crystalStructure) +{ + switch(crystalStructure) + { + case CrystalStructure::Hexagonal_High: + return 1.637833825000; + case CrystalStructure::Cubic_High: + return 1.096056815241; + case CrystalStructure::Hexagonal_Low: + return 3.141592653590; + case CrystalStructure::Cubic_Low: + return 1.570796326795; + case CrystalStructure::Triclinic: + return 3.141592653590; + case CrystalStructure::Monoclinic: + return 3.141592653590; + case CrystalStructure::OrthoRhombic: + return 2.094395102393; + case CrystalStructure::Tetragonal_Low: + return 3.141592653590; + case CrystalStructure::Tetragonal_High: + return 1.717771517458; + case CrystalStructure::Trigonal_Low: + return 3.141592653590; + case CrystalStructure::Trigonal_High: + return 1.823476581937; + default: + throw std::invalid_argument("random_angle_distribution::MaxMisorientationAngle: unknown crystal structure"); + } +} + +std::vector Compute(uint32_t crystalStructure, const std::vector& omega) +{ + const double maxAngle = MaxMisorientationAngle(crystalStructure); + + std::vector ad(omega.size(), 0.0); + + for(size_t i = 0; i < omega.size(); i++) + { + if(omega[i] > maxAngle) + { + ad[i] = 0.0; + continue; + } + + // omega/2 is in [0, maxAngle/2] with maxAngle <= pi, so mathematically + // rmag >= 0; std::abs guards the omega==maxAngle==pi boundary, where + // omega/2 sits right at tan()'s pole at pi/2 and floating-point rounding + // can push it a hair past the pole, flipping tan()'s sign. + const double rmag = std::abs(std::tan(omega[i] / 2.0)); + double xchi = 1.0; + + switch(crystalStructure) + { + case CrystalStructure::Hexagonal_Low: // 6/m + xchi = ChiCn(rmag, 6); + break; + case CrystalStructure::Tetragonal_Low: // 4/m + xchi = ChiCn(rmag, 4); + break; + case CrystalStructure::Trigonal_Low: // -3 + xchi = ChiCn(rmag, 3); + break; + case CrystalStructure::Monoclinic: // 2/m + xchi = ChiCn(rmag, 2); + break; + case CrystalStructure::Hexagonal_High: // 6/mmm + xchi = ChiDnh(rmag, 6); + break; + case CrystalStructure::Tetragonal_High: // 4/mmm + xchi = ChiDnh(rmag, 4); + break; + case CrystalStructure::Trigonal_High: // -3m + xchi = ChiDnh(rmag, 3); + break; + case CrystalStructure::OrthoRhombic: // mmm + xchi = ChiDnh(rmag, 2); + break; + case CrystalStructure::Cubic_Low: // m-3 + xchi = ChiM3(rmag); + break; + case CrystalStructure::Cubic_High: // m-3m + xchi = ChiM3m(rmag); + break; + case CrystalStructure::Triclinic: // -1 + xchi = 1.0; + break; + default: + throw std::invalid_argument("random_angle_distribution::Compute: unknown crystal structure"); + } + + // MTEX's `2 * numSym(cs)` prefactor is dropped here: it is a constant + // multiplier that cancels out in the unit-mean normalization below. + ad[i] = xchi * std::sin(omega[i] / 2.0) * std::sin(omega[i] / 2.0); + } + + double sum = 0.0; + size_t count = 0; + for(size_t i = 0; i < omega.size(); i++) + { + if(omega[i] <= maxAngle) + { + sum += ad[i]; + count++; + } + } + const double mean = (count > 0) ? (sum / static_cast(count)) : 0.0; + + for(size_t i = 0; i < ad.size(); i++) + { + if(omega[i] > maxAngle || mean == 0.0) + { + ad[i] = 0.0; + continue; + } + ad[i] /= mean; + if(ad[i] < 0.0) + { + ad[i] = 0.0; + } + } + + return ad; +} +} // namespace random_angle_distribution +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/RandomAngleDistribution.h b/Source/EbsdLib/Texture/RandomAngleDistribution.h new file mode 100644 index 00000000..3876d076 --- /dev/null +++ b/Source/EbsdLib/Texture/RandomAngleDistribution.h @@ -0,0 +1,61 @@ +/* ============================================================================ + * 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/EbsdLib.h" + +#include +#include + +namespace ebsdlib +{ +namespace random_angle_distribution +{ +/** + * @brief Maximum rotation angle (radians) of the fundamental region for the Laue group. + * Values generated from MTEX 6.1.0 fundamentalRegion(cs).maxAngle. + */ +EbsdLib_EXPORT double MaxMisorientationAngle(uint32_t crystalStructure); + +/** + * @brief Misorientation-angle distribution of the uniform (random) ODF. + * Port of MTEX geometry/@symmetry/calcAngleDistribution.m. Result is normalized + * to unit mean and zero-clamped. omega values beyond MaxMisorientationAngle get 0. + * Throws std::invalid_argument for UnknownCrystalStructure. + */ +EbsdLib_EXPORT std::vector Compute(uint32_t crystalStructure, const std::vector& omega); +} // namespace random_angle_distribution +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.cpp b/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.cpp new file mode 100644 index 00000000..0ecbd18b --- /dev/null +++ b/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.cpp @@ -0,0 +1,85 @@ +/* ============================================================================ + * 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 "SO3DeLaValleePoussinKernel.h" + +#include "EbsdLib/Math/EbsdLibMath.h" + +#include +#include + +namespace +{ +// std::beta is unavailable on libc++; use lgamma +double BetaFunction(double a, double b) +{ + return std::exp(std::lgamma(a) + std::lgamma(b) - std::lgamma(a + b)); +} +} // namespace + +namespace ebsdlib +{ +SO3DeLaValleePoussinKernel::SO3DeLaValleePoussinKernel(double halfwidthRadians) +: m_Halfwidth(halfwidthRadians) +{ + m_Kappa = 0.5 * std::log(0.5) / std::log(std::cos(halfwidthRadians / 2.0)); + m_C = BetaFunction(1.5, 0.5) / BetaFunction(1.5, m_Kappa + 0.5); +} + +double SO3DeLaValleePoussinKernel::kappa() const +{ + return m_Kappa; +} + +double SO3DeLaValleePoussinKernel::constant() const +{ + return m_C; +} + +double SO3DeLaValleePoussinKernel::halfwidth() const +{ + return m_Halfwidth; +} + +double SO3DeLaValleePoussinKernel::evaluate(double cosHalfOmega) const +{ + return m_C * std::pow(cosHalfOmega, 2.0 * m_Kappa); +} + +double SO3DeLaValleePoussinKernel::cutoffAngle() const +{ + return std::min(constants::k_PiD, 3.5 * m_Halfwidth); +} +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.h b/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.h new file mode 100644 index 00000000..ccb43a06 --- /dev/null +++ b/Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.h @@ -0,0 +1,76 @@ +/* ============================================================================ + * 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/EbsdLib.h" + +namespace ebsdlib +{ +/** + * @brief De la Vallee Poussin kernel on SO(3). Port of MTEX SO3DeLaValleePoussinKernel. + * + * K(omega) = C * cos(omega/2)^(2*kappa) with + * kappa = ln(0.5) / (2 * ln(cos(halfwidth/2))) + * C = B(1.5, 0.5) / B(1.5, kappa + 0.5) + * The kernel integrates to 1 over SO(3) with normalized Haar measure, so a + * weights-sum-to-one mixture of kernels is a normalized density (uniform == 1). + */ +class EbsdLib_EXPORT SO3DeLaValleePoussinKernel +{ +public: + explicit SO3DeLaValleePoussinKernel(double halfwidthRadians); + + double kappa() const; + double constant() const; + double halfwidth() const; + + /** + * @brief Evaluate the kernel. + * @param cosHalfOmega cos(omega/2); pass the absolute quaternion dot product. + */ + double evaluate(double cosHalfOmega) const; + + /** + * @brief Angle beyond which the kernel is treated as zero: min(pi, 3.5*halfwidth). + */ + double cutoffAngle() const; + +private: + double m_Halfwidth = 0.0; + double m_Kappa = 90.0; + double m_C = 0.0; +}; +} // namespace ebsdlib diff --git a/Source/EbsdLib/Texture/SourceList.cmake b/Source/EbsdLib/Texture/SourceList.cmake index 2bc936e0..01c871f7 100644 --- a/Source/EbsdLib/Texture/SourceList.cmake +++ b/Source/EbsdLib/Texture/SourceList.cmake @@ -4,10 +4,16 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TexturePreset.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/Texture.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/StatsGen.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/SO3DeLaValleePoussinKernel.h + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/RandomAngleDistribution.h + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/MisorientationKDE.h ) set(EbsdLib_${DIR_NAME}_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TexturePreset.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/SO3DeLaValleePoussinKernel.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/RandomAngleDistribution.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/MisorientationKDE.cpp ) #cmp_IDE_SOURCE_PROPERTIES("Common" "${EbsdLib_Texture_HDRS}" "${EbsdLib_Texture_SRCS}" "0") diff --git a/Source/EbsdLib/Utilities/inipp.h b/Source/EbsdLib/Utilities/inipp.h index 313a77ad..b1da437d 100644 --- a/Source/EbsdLib/Utilities/inipp.h +++ b/Source/EbsdLib/Utilities/inipp.h @@ -197,9 +197,9 @@ class Ini static const int max_interpolation_depth = 10; Ini() - : format(std::make_shared>()){}; + : format(std::make_shared>()) {}; Ini(std::shared_ptr> fmt) - : format(fmt){}; + : format(fmt) {}; void generate(std::basic_ostream& os) const { diff --git a/Source/EbsdLib/Utilities/wlenthe_orientation_coloring.hpp b/Source/EbsdLib/Utilities/wlenthe_orientation_coloring.hpp index 3edf2396..2b353748 100644 --- a/Source/EbsdLib/Utilities/wlenthe_orientation_coloring.hpp +++ b/Source/EbsdLib/Utilities/wlenthe_orientation_coloring.hpp @@ -215,7 +215,7 @@ void dihedralTriangle(T const* const n, T* const nTri) template bool cubicLowTriangle(T const* const n, T* const nTri) { - std::transform(n, n + 3, nTri, (T(*)(T)) & std::fabs); + std::transform(n, n + 3, nTri, (T(*)(T))&std::fabs); if(nTri[0] >= nTri[1]) { if(nTri[0] > nTri[2]) @@ -237,7 +237,7 @@ bool cubicLowTriangle(T const* const n, T* const nTri) template void cubicTriangle(T const* const n, T* const nTri) { - std::transform(n, n + 3, nTri, (T(*)(T)) & std::fabs); + std::transform(n, n + 3, nTri, (T(*)(T))&std::fabs); std::sort(nTri, nTri + 3); std::swap(nTri[0], nTri[1]); } diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index 7c0e0647..b5b5ef4f 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -40,6 +40,10 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/TexturePresetTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/PhaseTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/LaueOpsTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/MdfFZRodTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/DeLaValleePoussinKernelTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/RandomAngleDistributionTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/MisorientationKDETest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/ModifiedLambertProjection3DTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/ModifiedLambertProjectionTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/OrientationTransformationTest.cpp diff --git a/Source/Test/DeLaValleePoussinKernelTest.cpp b/Source/Test/DeLaValleePoussinKernelTest.cpp new file mode 100644 index 00000000..3471af4d --- /dev/null +++ b/Source/Test/DeLaValleePoussinKernelTest.cpp @@ -0,0 +1,49 @@ +#include + +#include "EbsdLib/Math/EbsdLibMath.h" +#include "EbsdLib/Texture/SO3DeLaValleePoussinKernel.h" + +#include + +using namespace ebsdlib; + +namespace +{ +constexpr double k_DegToRad = ebsdlib::constants::k_PiD / 180.0; +} + +// ----------------------------------------------------------------------------- +// SO3DeLaValleePoussinKernel is a direct port of MTEX 6.1.0's +// SO3DeLaValleePoussinKernel class (SO3Fun/SO3KernelFunctions/SO3DeLaValleePoussinKernel.m). +// The reference constants below were generated with: +// +// psi=SO3DeLaValleePoussinKernel('halfwidth',10*degree); +// fprintf('kappa10=%.12f C10=%.12f\n',psi.kappa,psi.C); +// psi5=SO3DeLaValleePoussinKernel('halfwidth',5*degree); +// fprintf('kappa5=%.12f C5=%.12f\n',psi5.kappa,psi5.C); +// fprintf('eval10_at20deg=%.12f\n', psi.eval(cos(20*degree/2))); +// +// MTEX's eval(psi, co2) takes co2 = cos(omega/2), matching evaluate(cosHalfOmega) +// here, so psi.eval(cos(20*degree/2)) is K(omega=20deg) for the halfwidth=10deg kernel. +TEST_CASE("SO3DeLaValleePoussinKernel matches MTEX", "[DeLaValleePoussinKernel]") +{ + SO3DeLaValleePoussinKernel psi10(10.0 * k_DegToRad); + REQUIRE(psi10.kappa() == Approx(90.903105993155).epsilon(1.0e-9)); + REQUIRE(psi10.constant() == Approx(1555.219446506688).epsilon(1.0e-9)); + REQUIRE(psi10.evaluate(std::cos(20.0 * k_DegToRad / 2.0)) == Approx(96.171329109334).epsilon(1.0e-9)); + + SO3DeLaValleePoussinKernel psi5(5.0 * k_DegToRad); + REQUIRE(psi5.kappa() == Approx(363.959328000956).epsilon(1.0e-9)); + REQUIRE(psi5.constant() == Approx(12345.110683440258).epsilon(1.0e-9)); + + // halfwidth round-trip: hw = 2*acos(0.5^(1/(2*kappa))) + REQUIRE(2.0 * std::acos(std::pow(0.5, 1.0 / (2.0 * psi10.kappa()))) == Approx(10.0 * k_DegToRad).epsilon(1.0e-12)); + + // halfwidth definition: K(hw) == K(0)/2 + REQUIRE(psi10.evaluate(std::cos(10.0 * k_DegToRad / 2.0)) == Approx(psi10.evaluate(1.0) / 2.0).epsilon(1.0e-9)); + + // cutoff = min(pi, 3.5*hw) + REQUIRE(psi10.cutoffAngle() == Approx(3.5 * 10.0 * k_DegToRad)); + SO3DeLaValleePoussinKernel psiWide(60.0 * k_DegToRad); + REQUIRE(psiWide.cutoffAngle() == Approx(ebsdlib::constants::k_PiD)); +} diff --git a/Source/Test/LaueOpsTest.cpp b/Source/Test/LaueOpsTest.cpp index 69103121..5c3d1898 100644 --- a/Source/Test/LaueOpsTest.cpp +++ b/Source/Test/LaueOpsTest.cpp @@ -819,3 +819,29 @@ TEST_CASE("ebsdlib::LaueOpsTest::IPFColor_SSTCorners", "[EbsdLib][LaueOpsTest]") checkDominant(ipfColorAtEtaChi(TrigonalLowOps(), -1.0 * oneDeg, chi90), 2, "TrigLow eta~0 (blue)"); } } + +// ----------------------------------------------------------------------------- +// The seven Laue classes that enumerate no slip systems used to set schmidfactor and +// slipsys but leave angleComps UNTOUCHED. That is silent: a caller which hoists one +// angleComps[2] buffer outside a per-Feature loop -- the natural way to write the loop -- +// got the PREVIOUS Feature's angle components attributed to a Feature of one of these +// classes, while schmidfactor correctly read 0. Pre-poisoning the buffer is what makes it +// observable: an overload that does not write angleComps leaves the poison in place. +TEST_CASE("ebsdlib::LaueOpsTest::SchmidFactorStubClassesDefineEveryOutput", "[EbsdLib][LaueOpsTest]") +{ + const std::vector stubOps = {TrigonalOps::New(), TrigonalLowOps::New(), TetragonalOps::New(), TetragonalLowOps::New(), + OrthoRhombicOps::New(), MonoclinicOps::New(), TriclinicOps::New()}; + for(const auto& op : stubOps) + { + double load[3] = {1.0, 2.0, 3.0}; + double schmidFactor = 7.0; + double angleComps[2] = {7.0, 9.0}; + int slipSystem = 7; + op->getSchmidFactorAndSS(load, schmidFactor, angleComps, slipSystem); + INFO(op->getNameOfClass()); + CHECK(schmidFactor == 0.0); + CHECK(slipSystem == 0); + CHECK(angleComps[0] == 0.0); + CHECK(angleComps[1] == 0.0); + } +} diff --git a/Source/Test/MdfFZRodTest.cpp b/Source/Test/MdfFZRodTest.cpp new file mode 100644 index 00000000..3b3d7cfd --- /dev/null +++ b/Source/Test/MdfFZRodTest.cpp @@ -0,0 +1,321 @@ +#include + +#include "EbsdLib/Core/EbsdLibConstants.h" +#include "EbsdLib/LaueOps/CubicLowOps.h" +#include "EbsdLib/LaueOps/CubicOps.h" +#include "EbsdLib/LaueOps/HexagonalLowOps.h" +#include "EbsdLib/LaueOps/HexagonalOps.h" +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/LaueOps/MonoclinicOps.h" +#include "EbsdLib/LaueOps/OrthoRhombicOps.h" +#include "EbsdLib/LaueOps/TetragonalLowOps.h" +#include "EbsdLib/LaueOps/TetragonalOps.h" +#include "EbsdLib/LaueOps/TriclinicOps.h" +#include "EbsdLib/LaueOps/TrigonalLowOps.h" +#include "EbsdLib/LaueOps/TrigonalOps.h" +#include "EbsdLib/Orientation/AxisAngle.hpp" +#include "EbsdLib/Orientation/Rodrigues.hpp" + +#include +#include +#include +#include +#include + +using namespace ebsdlib; + +// ----------------------------------------------------------------------------- +// getMDFFZRod() must map every symmetry-equivalent description of a misorientation +// to a single canonical Rodrigues vector in the MDF fundamental zone. Two +// misorientation axes n and n' are equivalent when +// +// n' = R(S_i) * n (conjugation by a rotational symmetry operator S_i), or +// n' = -R(S_i) * n (the above combined with switching symmetry: a boundary +// has no preferred direction, so a misorientation and its +// inverse -- which negates the axis -- are the same), and +// +// the misorientation angle w is invariant under both. These tests verify, for all +// eleven Laue classes: +// +// 1. Completeness: every equivalent axis folds to the same canonical rod. +// 2. Angle preservation: the rotation angle survives the fold unchanged. +// 3. Idempotence: the canonical rod is a fixed point of getMDFFZRod(). +// 4. Distinctness: axes that are NOT symmetry-equivalent fold to different rods +// (guards against over-folding; each pair below was wrongly merged by a +// historical implementation). +// +// The rotation angles used are kept below 30 degrees -- half the smallest nonzero +// symmetry rotation (60 degrees for the 6-fold) -- so the minimum-angle +// representative selection inside getMDFFZRod() is the identity for every variant +// and the tests isolate the axis-fold logic itself. +namespace +{ +std::array matVec(const Matrix3X3D& m, const std::array& v) +{ + return {m(0, 0) * v[0] + m(0, 1) * v[1] + m(0, 2) * v[2], m(1, 0) * v[0] + m(1, 1) * v[1] + m(1, 2) * v[2], m(2, 0) * v[0] + m(2, 1) * v[1] + m(2, 2) * v[2]}; +} + +std::array normalize(std::array v) +{ + const double mag = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + return {v[0] / mag, v[1] / mag, v[2] / mag}; +} + +RodriguesDType rodFromAxisAngle(const std::array& axis, double angle) +{ + return AxisAngleDType(axis[0], axis[1], axis[2], angle).toRodrigues(); +} + +// Compare two Rodrigues vectors through their axis-angle representation. +void checkSameRod(const RodriguesDType& computed, const RodriguesDType& expected, const std::string& what) +{ + AxisAngleDType axComputed = computed.toAxisAngle(); + AxisAngleDType axExpected = expected.toAxisAngle(); + INFO(what << " computed=(" << axComputed[0] << ", " << axComputed[1] << ", " << axComputed[2] << ", " << axComputed[3] << ")" + << " expected=(" << axExpected[0] << ", " << axExpected[1] << ", " << axExpected[2] << ", " << axExpected[3] << ")"); + CHECK(axComputed[0] == Approx(axExpected[0]).margin(1.0e-5)); + CHECK(axComputed[1] == Approx(axExpected[1]).margin(1.0e-5)); + CHECK(axComputed[2] == Approx(axExpected[2]).margin(1.0e-5)); + CHECK(axComputed[3] == Approx(axExpected[3]).margin(1.0e-6)); +} + +bool rodsDiffer(const RodriguesDType& a, const RodriguesDType& b) +{ + AxisAngleDType axA = a.toAxisAngle(); + AxisAngleDType axB = b.toAxisAngle(); + double dist = 0.0; + for(size_t i = 0; i < 3; i++) + { + dist += std::fabs(axA[i] - axB[i]); + } + return dist > 1.0e-3; +} + +std::array randomAxis(std::mt19937& gen) +{ + std::normal_distribution dist(0.0, 1.0); + while(true) + { + std::array v = {dist(gen), dist(gen), dist(gen)}; + const double mag = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + if(mag > 1.0e-3) + { + return {v[0] / mag, v[1] / mag, v[2] / mag}; + } + } +} + +// An axis at basal-plane azimuth (degrees) with the given z component. +std::array axisAtAzimuth(double azimuthDegrees, double z) +{ + const double azimuth = azimuthDegrees * ebsdlib::constants::k_PiOver180D; + const double mag = std::sqrt(1.0 - z * z); + return {mag * std::cos(azimuth), mag * std::sin(azimuth), z}; +} + +// Verify two inequivalent axes fold to different canonical rods. +void checkDistinct(const LaueOps& ops, const std::array& axisA, const std::array& axisB, double angle, const std::string& what) +{ + RodriguesDType foldA = ops.getMDFFZRod(rodFromAxisAngle(axisA, angle)); + RodriguesDType foldB = ops.getMDFFZRod(rodFromAxisAngle(axisB, angle)); + AxisAngleDType axA = foldA.toAxisAngle(); + AxisAngleDType axB = foldB.toAxisAngle(); + INFO(what << " foldA=(" << axA[0] << ", " << axA[1] << ", " << axA[2] << ")" + << " foldB=(" << axB[0] << ", " << axB[1] << ", " << axB[2] << ")"); + CHECK(rodsDiffer(foldA, foldB)); +} +} // namespace + +// ----------------------------------------------------------------------------- +TEST_CASE("ebsdlib::MdfFZRodTest::EquivalentAxesCollapse", "[EbsdLib][MdfFZRodTest]") +{ + constexpr size_t k_NumRandomAxes = 25; + const std::vector k_Angles = {2.0 * ebsdlib::constants::k_PiOver180D, 25.0 * ebsdlib::constants::k_PiOver180D}; + + auto allOps = LaueOps::GetAllOrientationOps(); + size_t classesChecked = 0; + + for(size_t laueIdx = 0; laueIdx < 11; laueIdx++) + { + const LaueOps& ops = *allOps[laueIdx]; + const std::string className = ops.getNameOfClass(); + const size_t numSym = ops.getNumSymOps(); + + SECTION(className) + { + std::mt19937 gen(20260722); + for(double angle : k_Angles) + { + for(size_t trial = 0; trial < k_NumRandomAxes; trial++) + { + const std::array axis = randomAxis(gen); + const RodriguesDType reference = ops.getMDFFZRod(rodFromAxisAngle(axis, angle)); + + // Angle preservation: the fold must not alter the misorientation angle + { + AxisAngleDType axRef = reference.toAxisAngle(); + INFO(className << " angle preservation, trial " << trial); + CHECK(axRef[3] == Approx(angle).margin(1.0e-6)); + } + + // Idempotence: the canonical rod must be a fixed point + checkSameRod(ops.getMDFFZRod(reference), reference, className + " idempotence"); + + // Completeness: every conjugated and/or switched axis folds to the reference + for(size_t symIdx = 0; symIdx < numSym; symIdx++) + { + const std::array conjugated = matVec(ops.getMatSymOpD(symIdx), axis); + const std::array switched = {-conjugated[0], -conjugated[1], -conjugated[2]}; + + checkSameRod(ops.getMDFFZRod(rodFromAxisAngle(conjugated, angle)), reference, className + " conjugation by op " + std::to_string(symIdx)); + checkSameRod(ops.getMDFFZRod(rodFromAxisAngle(switched, angle)), reference, className + " switching + op " + std::to_string(symIdx)); + } + } + } + classesChecked++; + } + } +} + +// ----------------------------------------------------------------------------- +// Axes lying exactly on symmetry elements (the c-axis, a 2-fold axis, a body +// diagonal) sit on fold-sector boundaries. The fold must handle them without +// producing NaN and must still collapse their +/- variants. +TEST_CASE("ebsdlib::MdfFZRodTest::BoundaryAxes", "[EbsdLib][MdfFZRodTest]") +{ + const double angle = 10.0 * ebsdlib::constants::k_PiOver180D; + const double k_Root3Inverse = 1.0 / std::sqrt(3.0); + const std::vector> boundaryAxes = { + {0.0, 0.0, 1.0}, + {1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {k_Root3Inverse, k_Root3Inverse, k_Root3Inverse}, + }; + + auto allOps = LaueOps::GetAllOrientationOps(); + + for(size_t laueIdx = 0; laueIdx < 11; laueIdx++) + { + const LaueOps& ops = *allOps[laueIdx]; + const std::string className = ops.getNameOfClass(); + + SECTION(className) + { + for(const std::array& axis : boundaryAxes) + { + const RodriguesDType folded = ops.getMDFFZRod(rodFromAxisAngle(axis, angle)); + AxisAngleDType axFolded = folded.toAxisAngle(); + + INFO(className << " boundary axis (" << axis[0] << ", " << axis[1] << ", " << axis[2] << ")"); + for(size_t i = 0; i < 4; i++) + { + CHECK(std::isfinite(axFolded[i])); + } + CHECK(axFolded[3] == Approx(angle).margin(1.0e-6)); + + // The negated axis is the switched (inverse) misorientation: must collapse + const std::array negated = {-axis[0], -axis[1], -axis[2]}; + checkSameRod(ops.getMDFFZRod(rodFromAxisAngle(negated, angle)), folded, className + " negated boundary axis"); + } + } + } +} + +// ----------------------------------------------------------------------------- +// Axes that are NOT symmetry-equivalent must fold to DIFFERENT canonical rods. +// Every pair below was wrongly merged by a historical implementation of +// getMDFFZRod(), so these are regression guards against over-folding. +TEST_CASE("ebsdlib::MdfFZRodTest::InequivalentAxesStayDistinct", "[EbsdLib][MdfFZRodTest]") +{ + const double angle = 10.0 * ebsdlib::constants::k_PiOver180D; + + SECTION("CubicLowOps: m-3 has no 2-component transposition") + { + // The tetrahedral group's <111> 3-folds give only cyclic permutations of the + // axis components; a transposition of the two smaller components is not a + // symmetry (the old full descending sort merged these). + CubicLowOps ops; + checkDistinct(ops, normalize({0.8, 0.5, 0.33}), normalize({0.8, 0.33, 0.5}), angle, "CubicLow transposition"); + } + + SECTION("TetragonalOps: 422 has no 3-fold") + { + TetragonalOps ops; + checkDistinct(ops, normalize({0.8, 0.5, 0.33}), normalize({0.33, 0.8, 0.5}), angle, "Tetragonal cyclic permutation"); + } + + SECTION("TetragonalLowOps: 4/m has no in-plane mirror") + { + // Azimuth +20 and -20 degrees are mirror images; 4/m has no in-plane 2-folds + // so they are inequivalent (the old fabs fold merged them). + TetragonalLowOps ops; + checkDistinct(ops, axisAtAzimuth(20.0, 0.7), axisAtAzimuth(-20.0, 0.7), angle, "TetragonalLow azimuth mirror"); + } + + SECTION("HexagonalLowOps: 6/m has no in-plane mirror") + { + HexagonalLowOps ops; + checkDistinct(ops, axisAtAzimuth(20.0, 0.7), axisAtAzimuth(-20.0, 0.7), angle, "HexagonalLow azimuth mirror"); + } + + SECTION("TrigonalLowOps: -3 has no in-plane mirror") + { + TrigonalLowOps ops; + checkDistinct(ops, axisAtAzimuth(20.0, 0.7), axisAtAzimuth(-20.0, 0.7), angle, "TrigonalLow azimuth mirror"); + } + + SECTION("TrigonalOps: mirror lines sit at 30+60k degrees, not 60k") + { + // Azimuths 70 and 50 degrees are NOT related by the 32 group's mirrors at + // 30/90/150 degrees (the old fan-fold with boundaries at 60k merged them). + TrigonalOps ops; + checkDistinct(ops, axisAtAzimuth(70.0, 0.7), axisAtAzimuth(50.0, 0.7), angle, "Trigonal 70 vs 50 degrees"); + } + + SECTION("TriclinicOps: only switching symmetry") + { + TriclinicOps ops; + checkDistinct(ops, normalize({0.5, 0.6, 0.62}), normalize({-0.5, 0.6, 0.62}), angle, "Triclinic single-component flip"); + } + + SECTION("MonoclinicOps: single-component flip of n1 is not a symmetry") + { + MonoclinicOps ops; + checkDistinct(ops, normalize({0.5, 0.6, 0.62}), normalize({-0.5, 0.6, 0.62}), angle, "Monoclinic single-component flip"); + } + + SECTION("OrthoRhombicOps: mmm has no permutation symmetry") + { + OrthoRhombicOps ops; + checkDistinct(ops, normalize({0.8, 0.5, 0.33}), normalize({0.5, 0.8, 0.33}), angle, "OrthoRhombic transposition"); + } +} + +// ----------------------------------------------------------------------------- +// The misorientation angle discarded by historical CubicLowOps / TrigonalOps / +// TrigonalLowOps implementations: verify explicitly for every class that a +// nonzero input angle never comes back as zero. +TEST_CASE("ebsdlib::MdfFZRodTest::AngleNeverDiscarded", "[EbsdLib][MdfFZRodTest]") +{ + const double angle = 15.0 * ebsdlib::constants::k_PiOver180D; + + auto allOps = LaueOps::GetAllOrientationOps(); + + for(size_t laueIdx = 0; laueIdx < 11; laueIdx++) + { + const LaueOps& ops = *allOps[laueIdx]; + const std::string className = ops.getNameOfClass(); + + SECTION(className) + { + std::mt19937 gen(42); + for(size_t trial = 0; trial < 5; trial++) + { + const std::array axis = randomAxis(gen); + AxisAngleDType axFolded = ops.getMDFFZRod(rodFromAxisAngle(axis, angle)).toAxisAngle(); + INFO(className << " trial " << trial); + CHECK(axFolded[3] == Approx(angle).margin(1.0e-6)); + } + } + } +} diff --git a/Source/Test/MisorientationKDETest.cpp b/Source/Test/MisorientationKDETest.cpp new file mode 100644 index 00000000..122dd0e0 --- /dev/null +++ b/Source/Test/MisorientationKDETest.cpp @@ -0,0 +1,405 @@ +#include + +#include "EbsdLib/Core/EbsdLibConstants.h" +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/Math/EbsdLibMath.h" +#include "EbsdLib/Orientation/AxisAngle.hpp" +#include "EbsdLib/Orientation/Quaternion.hpp" +#include "EbsdLib/Orientation/Rodrigues.hpp" +#include "EbsdLib/Texture/MisorientationKDE.h" +#include "EbsdLib/Texture/RandomAngleDistribution.h" +#include "EbsdLib/Texture/SO3DeLaValleePoussinKernel.h" + +#include +#include +#include +#include + +using namespace ebsdlib; + +namespace +{ +constexpr double k_DegToRad = ebsdlib::constants::k_PiOver180D; + +// Build a unit quaternion from a (not necessarily unit) axis and an angle in radians. +QuatD quatFromAxisAngle(double ax, double ay, double az, double angleRadians) +{ + const double mag = std::sqrt(ax * ax + ay * ay + az * az); + const double nx = ax / mag; + const double ny = ay / mag; + const double nz = az / mag; + const double s = std::sin(angleRadians / 2.0); + return QuatD(nx * s, ny * s, nz * s, std::cos(angleRadians / 2.0)); +} + +} // namespace + +// ----------------------------------------------------------------------------- +// Triclinic has a single (identity) symmetry operator, so the crystal-symmetry +// average is trivial and the kernel density around one isolated misorientation +// reduces to the raw De la Vallee Poussin kernel scaled by the 0.5 antipodal +// normalization factor. Its modal height is K(0) / 2 = psi.evaluate(1.0) / 2, it +// falls to a quarter peak a halfwidth away, and it is exactly zero beyond the cutoff. +// (The 0.5 factor was pinned by the MTEX cross-check in CubicAngleCurveVsMTEX.) +TEST_CASE("ebsdlib::MisorientationKDE::SingleCenterTriclinic", "[EbsdLib][MisorientationKDE]") +{ + auto opsList = ebsdlib::LaueOps::GetAllOrientationOps(); + auto ops = opsList[ebsdlib::CrystalStructure::Triclinic]; + const double hw = 10.0 * k_DegToRad; + ebsdlib::MisorientationKDE kde(ops, ebsdlib::CrystalStructure::Triclinic, hw); + + // center: 30 degrees about z + ebsdlib::QuatD c(0.0, 0.0, std::sin(15.0 * k_DegToRad), std::cos(15.0 * k_DegToRad)); + kde.addMisorientation(c, 3.0); // non-unit weight; must normalize to 1 + kde.finalize(); + REQUIRE(kde.totalWeight() == Approx(3.0)); + + ebsdlib::SO3DeLaValleePoussinKernel psi(hw); + + // The KDE center is the bin's weighted circular-mean misorientation, which for a single + // observation is the (FZ-folded) input misorientation itself, so the modal peak sits at that + // misorientation rather than at the geometric bin center. + ebsdlib::QuatD center = ops->getMDFFZRod(c.toRodrigues()).toQuaternion(); + + // Modal peak: density at the center equals K(0) / 2 (the 0.5 antipodal factor). + REQUIRE(kde.evaluate(center) == Approx(psi.evaluate(1.0) / 2.0).epsilon(0.01)); + + // Rotate the center by hw about an orthogonal axis -> quarter peak (half of K(0)/2). + ebsdlib::QuatD dq(std::sin(hw / 2.0), 0.0, 0.0, std::cos(hw / 2.0)); + REQUIRE(kde.evaluate(dq * center) == Approx(psi.evaluate(1.0) / 4.0).epsilon(0.02)); + + // Beyond the cutoff -> exactly zero. + ebsdlib::QuatD far(0.0, std::sin(60.0 * k_DegToRad), 0.0, std::cos(60.0 * k_DegToRad)); + REQUIRE(kde.evaluate(far * center) == 0.0); +} + +// ----------------------------------------------------------------------------- +// Cubic (m-3m) exercises the full crystal-symmetry averaging and the antipodal +// (grain-exchange) folding: the density must be invariant under s_1 * q * s_2 for +// every symmetry pair, invariant under grain exchange q -> q^-1, and normalized so +// that its mean over the fundamental zone is close to 1. +TEST_CASE("ebsdlib::MisorientationKDE::CubicInvarianceAndMean", "[EbsdLib][MisorientationKDE]") +{ + auto ops = ebsdlib::LaueOps::GetAllOrientationOps()[ebsdlib::CrystalStructure::Cubic_High]; + ebsdlib::MisorientationKDE kde(ops, ebsdlib::CrystalStructure::Cubic_High, 10.0 * k_DegToRad); + + // Three arbitrary misorientations with unequal weights. + kde.addMisorientation(quatFromAxisAngle(0.0, 0.0, 1.0, 25.0 * k_DegToRad), 1.0); + kde.addMisorientation(quatFromAxisAngle(1.0, 1.0, 1.0, 40.0 * k_DegToRad), 2.0); + kde.addMisorientation(quatFromAxisAngle(0.0, 1.0, 2.0, 55.0 * k_DegToRad), 3.0); + kde.finalize(); + + REQUIRE(kde.totalWeight() == Approx(6.0)); + + // A generic query misorientation. + ebsdlib::QuatD query = quatFromAxisAngle(1.0, 2.0, 3.0, 33.0 * k_DegToRad); + const double reference = kde.evaluate(query); + REQUIRE(reference > 0.0); + + const size_t numSymOps = ops->getNumSymOps(); + + SECTION("crystal symmetry invariance: f(s1 * q * s2) == f(q)") + { + const std::array leftIdx = {0, 1, 5, 11}; + const std::array rightIdx = {0, 2, 7, 13}; + for(size_t li : leftIdx) + { + for(size_t ri : rightIdx) + { + if(li >= numSymOps || ri >= numSymOps) + { + continue; + } + ebsdlib::QuatD s1 = ops->getQuatSymOp(li); + ebsdlib::QuatD s2 = ops->getQuatSymOp(ri); + ebsdlib::QuatD equivalent = s1 * query * s2; + INFO("left op " << li << " right op " << ri); + CHECK(kde.evaluate(equivalent) == Approx(reference).epsilon(1.0e-6)); + } + } + } + + SECTION("grain-exchange invariance: f(q^-1) == f(q)") + { + CHECK(kde.evaluate(query.conjugate()) == Approx(reference).epsilon(1.0e-6)); + } + + SECTION("normalization: mean density over non-identity bin centers is O(1)") + { + std::vector densities = kde.evaluateAtBinCenters(); + REQUIRE(densities.size() == ops->getMDFSize()); + + double sum = 0.0; + size_t count = 0; + for(size_t binIndex = 0; binIndex < densities.size(); binIndex++) + { + ebsdlib::QuatD center = kde.binCenter(static_cast(binIndex)); + // Disorientation angle of the bin center from identity: omega = 2*acos(|w|). + const double cosHalf = std::fabs(center.w()); + const double omega = 2.0 * std::acos(std::min(1.0, cosHalf)); + if(omega > 1.0e-3) + { + sum += densities[binIndex]; + count++; + } + } + REQUIRE(count > 0); + const double mean = sum / static_cast(count); + INFO("mean density over " << count << " bins = " << mean); + CHECK(mean > 0.5); + CHECK(mean < 2.0); + } +} + +// ----------------------------------------------------------------------------- +// Numerical cross-check against MTEX 6.1.0. The same 3-center cubic KDE +// (weights 1,2,3) is built, and computeAngleCurve(200) is compared against +// MTEX's calcDensity(...,'exact') -> calcAngleDistribution reference at 20 +// sampled angles. This is the authority that pins the KDE's absolute scale. +// +// MTEX reference generated with: +// cs=crystalSymmetry('m-3m'); +// ax=[vector3d(0,0,1), vector3d(1,1,1)/norm(vector3d(1,1,1)), vector3d(0,1,2)/norm(vector3d(0,1,2))]; +// om=[25 40 55]*degree; mori=orientation('axis',ax,'angle',om,cs,cs); w=[1 2 3]; +// mdf=calcDensity(mori,'weights',w,'halfwidth',10*degree,'exact'); +// [d,omega]=calcAngleDistribution(mdf); +// MTEX omega(k) == maxAngle*(k-1)/199, so MTEX index k maps to C++ curve index k-1. +TEST_CASE("ebsdlib::MisorientationKDE::CubicAngleCurveVsMTEX", "[EbsdLib][MisorientationKDE]") +{ + auto ops = ebsdlib::LaueOps::GetAllOrientationOps()[ebsdlib::CrystalStructure::Cubic_High]; + ebsdlib::MisorientationKDE kde(ops, ebsdlib::CrystalStructure::Cubic_High, 10.0 * k_DegToRad); + kde.addMisorientation(quatFromAxisAngle(0.0, 0.0, 1.0, 25.0 * k_DegToRad), 1.0); + kde.addMisorientation(quatFromAxisAngle(1.0, 1.0, 1.0, 40.0 * k_DegToRad), 2.0); + kde.addMisorientation(quatFromAxisAngle(0.0, 1.0, 2.0, 55.0 * k_DegToRad), 3.0); + kde.finalize(); + + const size_t numPoints = 200; + ebsdlib::MisorientationKDE::AngleCurve curve = kde.computeAngleCurve(numPoints); + + REQUIRE(curve.Angles.size() == numPoints); + REQUIRE(curve.Density.size() == numPoints); + REQUIRE(curve.RandomDensity.size() == numPoints); + + // Angle grid endpoints. + const double maxAngle = ebsdlib::random_angle_distribution::MaxMisorientationAngle(ebsdlib::CrystalStructure::Cubic_High); + REQUIRE(curve.Angles.front() == Approx(0.0).margin(1.0e-12)); + REQUIRE(curve.Angles.back() == Approx(maxAngle)); + + // RandomDensity must match the analytic reference exactly. + std::vector expectedRandom = ebsdlib::random_angle_distribution::Compute(ebsdlib::CrystalStructure::Cubic_High, curve.Angles); + REQUIRE(expectedRandom.size() == numPoints); + for(size_t i = 0; i < numPoints; i++) + { + INFO("RandomDensity mismatch at index " << i); + CHECK(curve.RandomDensity[i] == Approx(expectedRandom[i]).margin(1.0e-12)); + } + + // 20 sampled MTEX (1-based index k, omega, density) reference pairs. + const std::array, 20> mtexRef = { + {{{1, 0.0000000000, 0.0000000000}}, {{11, 0.0550782319, 0.0020168845}}, {{21, 0.1101564638, 0.0111024402}}, {{31, 0.1652346958, 0.0366089735}}, {{41, 0.2203129277, 0.0928770466}}, + {{51, 0.2753911596, 0.1866315579}}, {{61, 0.3304693915, 0.3349090171}}, {{71, 0.3855476235, 0.5090465258}}, {{81, 0.4406258554, 0.7013347724}}, {{91, 0.4957040873, 0.9018987543}}, + {{101, 0.5507823192, 1.1704777051}}, {{111, 0.6058605511, 1.5494982458}}, {{121, 0.6609387831, 2.0076201248}}, {{131, 0.7160170150, 2.4918589687}}, {{141, 0.7710952469, 2.8632147418}}, + {{151, 0.8261734788, 2.9262518556}}, {{161, 0.8812517107, 2.1954465862}}, {{171, 0.9363299427, 1.1739709599}}, {{181, 0.9914081746, 0.5369126565}}, {{191, 1.0464864065, 0.2176428572}}}}; + + // Tolerance: the KDE math itself matches MTEX exactly. Evaluating the density at the + // *exact* (un-gridified) misorientation centers reproduces MTEX's mdf to a constant + // ratio of 2.000 at every angle (and mean(mdf) == 1 over SO(3)), which is what pinned + // the 0.5 antipodal normalization factor in MisorientationKDE::evaluate(). The residual + // pointwise deviation seen here (up to ~45% on the low-omega tail, ~18% on the rising + // flank) is entirely the Task 3 MDF-bin gridify: addMisorientation() snaps each center + // to its ~5-degree MDF fundamental-zone bin center, shifting the 25/40/55-degree inputs + // by 1-2.5 degrees, which redistributes the steep 10-degree-halfwidth kernel across + // omega. Removing the snap collapses the deviation to <0.1%. The tolerance below is + // therefore set to the gridify band (20% relative + 0.10 absolute floor), not to the + // KDE's intrinsic accuracy. See task-4-report.md for the exact-center evidence. + for(const std::array& ref : mtexRef) + { + const size_t idx = static_cast(std::lround(ref[0])) - 1; // 1-based -> 0-based + const double mtexOmega = ref[1]; + const double mtexDensity = ref[2]; + INFO("MTEX k=" << ref[0] << " omega=" << mtexOmega << " -> curve index " << idx << " angle=" << curve.Angles[idx] << " density=" << curve.Density[idx]); + CHECK(curve.Angles[idx] == Approx(mtexOmega).margin(1.0e-6)); + CHECK(curve.Density[idx] == Approx(mtexDensity).epsilon(0.20).margin(0.10)); + } +} + +// ----------------------------------------------------------------------------- +// Numerical cross-check against MTEX 6.1.0 for HEXAGONAL (6/mmm). This is the +// hexagonal analogue of CubicAngleCurveVsMTEX and closes the gap where the hex +// angle-curve extraction had never been pinned to MTEX (only the Mackenzie random +// reference had been). A 3-center hex KDE (weights 1,2,3) is built and +// computeAngleCurve(200) is compared against MTEX's calcDensity(...,'exact') -> +// calcAngleDistribution reference at 20 sampled angles. It also pins the hex +// MaxMisorientationAngle and the random-density reference against MTEX. +// +// MTEX reference generated with: +// cs=crystalSymmetry('6/mmm'); +// ax=[vector3d(0,0,1), vector3d(1,0,0), vector3d(1,1,1)/norm(vector3d(1,1,1))]; +// om=[20 50 80]*degree; mori=orientation('axis',ax,'angle',om,cs,cs); w=[1 2 3]; +// mdf=calcDensity(mori,'weights',w,'halfwidth',10*degree,'exact'); +// [d,omega]=calcAngleDistribution(mdf); +// MTEX omega(k) == maxAngle*(k-1)/199, so MTEX index k maps to C++ curve index k-1. +// The axes are cartesian vector3d in the crystal frame (x=a1, z=c), matching how +// the NX KDE builds the misorientation quaternion directly from a cartesian axis. +TEST_CASE("ebsdlib::MisorientationKDE::HexagonalAngleCurveVsMTEX", "[EbsdLib][MisorientationKDE]") +{ + auto ops = ebsdlib::LaueOps::GetAllOrientationOps()[ebsdlib::CrystalStructure::Hexagonal_High]; + ebsdlib::MisorientationKDE kde(ops, ebsdlib::CrystalStructure::Hexagonal_High, 10.0 * k_DegToRad); + const double invSqrt3 = 1.0 / std::sqrt(3.0); + kde.addMisorientation(quatFromAxisAngle(0.0, 0.0, 1.0, 20.0 * k_DegToRad), 1.0); + kde.addMisorientation(quatFromAxisAngle(1.0, 0.0, 0.0, 50.0 * k_DegToRad), 2.0); + kde.addMisorientation(quatFromAxisAngle(invSqrt3, invSqrt3, invSqrt3, 80.0 * k_DegToRad), 3.0); + kde.finalize(); + + const size_t numPoints = 200; + ebsdlib::MisorientationKDE::AngleCurve curve = kde.computeAngleCurve(numPoints); + + REQUIRE(curve.Angles.size() == numPoints); + REQUIRE(curve.Density.size() == numPoints); + REQUIRE(curve.RandomDensity.size() == numPoints); + + // Angle grid endpoints. The hex 6/mmm maximum misorientation angle is 93.84 degrees, + // which matches MTEX fundamentalRegion('6/mmm','6/mmm').maxAngle == 1.637833825 rad. + const double maxAngle = ebsdlib::random_angle_distribution::MaxMisorientationAngle(ebsdlib::CrystalStructure::Hexagonal_High); + REQUIRE(maxAngle == Approx(1.637833825).margin(1.0e-9)); + REQUIRE(curve.Angles.front() == Approx(0.0).margin(1.0e-12)); + REQUIRE(curve.Angles.back() == Approx(maxAngle)); + + // RandomDensity must match the analytic (Mackenzie) reference exactly; that reference + // was independently cross-checked against MTEX calcAngleDistribution('6/mmm') to 1e-12. + std::vector expectedRandom = ebsdlib::random_angle_distribution::Compute(ebsdlib::CrystalStructure::Hexagonal_High, curve.Angles); + REQUIRE(expectedRandom.size() == numPoints); + for(size_t i = 0; i < numPoints; i++) + { + INFO("RandomDensity mismatch at index " << i); + CHECK(curve.RandomDensity[i] == Approx(expectedRandom[i]).margin(1.0e-12)); + } + + // 20 sampled MTEX (1-based index k, omega, density) reference pairs. + const std::array, 20> mtexRef = { + {{{1, 0.0000000000, 0.0000000000}}, {{11, 0.0823032073, 0.0352621130}}, {{21, 0.1646064146, 0.1636836591}}, {{31, 0.2469096219, 0.4222072865}}, {{41, 0.3292128291, 0.7709140872}}, + {{51, 0.4115160364, 0.8567685815}}, {{61, 0.4938192437, 0.9308375971}}, {{71, 0.5761224510, 0.3756435243}}, {{81, 0.6584256583, 0.4738874020}}, {{91, 0.7407288656, 0.9702940778}}, + {{101, 0.8230320729, 1.5630454079}}, {{111, 0.9053352802, 2.2704391764}}, {{121, 0.9876384874, 2.7324230182}}, {{131, 1.0699416947, 3.0542566503}}, {{141, 1.1522449020, 2.6611603119}}, + {{151, 1.2345481093, 1.8815474011}}, {{161, 1.3168513166, 1.0179783863}}, {{171, 1.3991545239, 0.4290641078}}, {{181, 1.4814577312, 0.1270194532}}, {{191, 1.5637609384, 0.0338470058}}}}; + + // Tolerance mirrors CubicAngleCurveVsMTEX. Because the KDE centers each bin on the + // weighted circular mean of the misorientations that fell in it (not the geometric + // ~5-degree bin center), these three isolated hex misorientations are represented at + // essentially their exact positions, so the curve tracks MTEX's 'exact' reference + // closely across the whole angle range -- including the low-to-mid-angle band (14/19/38 + // degrees) that the earlier geometric-bin-center snap under-estimated by 23-32%. The + // band below (20% relative + 0.10 absolute floor) covers the residual axis-grid + // sampling noise and the modal peak angle matches MTEX exactly (curve index 131, + // ~61.3 degrees). + for(const std::array& ref : mtexRef) + { + const size_t idx = static_cast(std::lround(ref[0])) - 1; // 1-based -> 0-based + const double mtexOmega = ref[1]; + const double mtexDensity = ref[2]; + INFO("MTEX k=" << ref[0] << " omega=" << mtexOmega << " -> curve index " << idx << " angle=" << curve.Angles[idx] << " density=" << curve.Density[idx]); + CHECK(curve.Angles[idx] == Approx(mtexOmega).margin(1.0e-6)); + CHECK(curve.Density[idx] == Approx(mtexDensity).epsilon(0.20).margin(0.10)); + } + + // The modal peak of the measured hex MDF curve sits at ~61 degrees (MTEX index 131), + // well away from the hex random-reference (Mackenzie) maximum, confirming the curve + // carries the correlated-misorientation signal rather than collapsing to the reference. + size_t curveArgMax = 0; + for(size_t i = 1; i < curve.Density.size(); i++) + { + if(curve.Density[i] > curve.Density[curveArgMax]) + { + curveArgMax = i; + } + } + const double curvePeakDeg = curve.Angles[curveArgMax] / k_DegToRad; + INFO("hex angle-curve peak (deg) = " << curvePeakDeg); + CHECK(curvePeakDeg > 58.0); + CHECK(curvePeakDeg < 64.0); +} + +// ----------------------------------------------------------------------------- +// Regression guard for a correlated MDF whose true peak is at a KNOWN non-45-degree +// angle. A tight cluster of Sigma3 (60 degree / <111>) misorientations on top of a +// random background must produce an MDF whose bin-array peak folds to a 60/<111> +// misorientation AND whose angle-distribution curve peaks near 60 degrees -- clearly +// distinct from the ~45 degree cubic Mackenzie (random-reference) maximum. This is the +// discriminating case the earlier 45-degree bicrystal cross-check could not catch: a +// density that collapsed to the random distribution would still peak at ~45 and pass a +// weaker test, but fails here. +TEST_CASE("ebsdlib::MisorientationKDE::CorrelatedTwinPeaksAtSixty", "[EbsdLib][MisorientationKDE]") +{ + auto ops = ebsdlib::LaueOps::GetAllOrientationOps()[ebsdlib::CrystalStructure::Cubic_High]; + const double hw = 10.0 * k_DegToRad; + ebsdlib::MisorientationKDE kde(ops, ebsdlib::CrystalStructure::Cubic_High, hw); + + // Tight cluster of 60 degree / <111> Sigma3 twins. + const size_t numTwins = 600; + for(size_t i = 0; i < numTwins; i++) + { + kde.addMisorientation(quatFromAxisAngle(1.0, 1.0, 1.0, 60.0 * k_DegToRad), 1.0); + } + // Uniform-ish random background (deterministic) via Shoemake's method. + std::mt19937 gen(12345); + std::uniform_real_distribution uni(0.0, 1.0); + const size_t numRandom = 400; + for(size_t i = 0; i < numRandom; i++) + { + const double u1 = uni(gen); + const double u2 = uni(gen); + const double u3 = uni(gen); + const double s1 = std::sqrt(1.0 - u1); + const double s2 = std::sqrt(u1); + ebsdlib::QuatD q(s1 * std::sin(2.0 * constants::k_PiD * u2), s1 * std::cos(2.0 * constants::k_PiD * u2), s2 * std::sin(2.0 * constants::k_PiD * u3), s2 * std::cos(2.0 * constants::k_PiD * u3)); + kde.addMisorientation(q, 1.0); + } + kde.finalize(); + + // 1. The MDF bin-array peak folds to a 60 degree / <111> misorientation. + std::vector densities = kde.evaluateAtBinCenters(); + size_t argMax = 0; + for(size_t i = 1; i < densities.size(); i++) + { + if(densities[i] > densities[argMax]) + { + argMax = i; + } + } + double seed[3] = {0.5, 0.5, 0.5}; + ebsdlib::RodriguesDType peakRod = ops->determineRodriguesVector(seed, static_cast(argMax)); + ebsdlib::AxisAngleDType peakAxisAngle = peakRod.toAxisAngle(); + const double peakAngleDeg = peakAxisAngle[3] / k_DegToRad; + INFO("MDF peak angle (deg) = " << peakAngleDeg << " axis (" << peakAxisAngle[0] << ", " << peakAxisAngle[1] << ", " << peakAxisAngle[2] << ")"); + CHECK(peakAngleDeg > 56.0); + CHECK(peakAngleDeg < 63.0); + const double invSqrt3 = 1.0 / std::sqrt(3.0); + CHECK(std::fabs(std::fabs(peakAxisAngle[0]) - invSqrt3) < 0.1); + CHECK(std::fabs(std::fabs(peakAxisAngle[1]) - invSqrt3) < 0.1); + CHECK(std::fabs(std::fabs(peakAxisAngle[2]) - invSqrt3) < 0.1); + + // 2. The angle-distribution curve peaks near 60 degrees, not at the ~45 degree random peak. + ebsdlib::MisorientationKDE::AngleCurve curve = kde.computeAngleCurve(200); + size_t curveArgMax = 0; + for(size_t i = 1; i < curve.Density.size(); i++) + { + if(curve.Density[i] > curve.Density[curveArgMax]) + { + curveArgMax = i; + } + } + size_t randomArgMax = 0; + for(size_t i = 1; i < curve.RandomDensity.size(); i++) + { + if(curve.RandomDensity[i] > curve.RandomDensity[randomArgMax]) + { + randomArgMax = i; + } + } + const double curvePeakDeg = curve.Angles[curveArgMax] / k_DegToRad; + const double randomPeakDeg = curve.Angles[randomArgMax] / k_DegToRad; + INFO("angle-curve peak (deg) = " << curvePeakDeg << ", random-reference peak (deg) = " << randomPeakDeg); + CHECK(curvePeakDeg > 53.0); + CHECK(curvePeakDeg < 63.0); + // The random reference itself peaks near 45 degrees; the measured curve must be well above it. + CHECK(curvePeakDeg - randomPeakDeg > 5.0); +} diff --git a/Source/Test/RandomAngleDistributionTest.cpp b/Source/Test/RandomAngleDistributionTest.cpp new file mode 100644 index 00000000..bfdc7aeb --- /dev/null +++ b/Source/Test/RandomAngleDistributionTest.cpp @@ -0,0 +1,150 @@ +#include + +#include "EbsdLib/Core/EbsdLibConstants.h" +#include "EbsdLib/Texture/RandomAngleDistribution.h" + +#include +#include + +using namespace ebsdlib; + +namespace +{ +// ----------------------------------------------------------------------------- +// RandomAngleDistribution is a direct port of MTEX 6.1.0's +// geometry/@symmetry/calcAngleDistribution.m (lines 44-218). The reference +// constants below were generated with: +// +// names={'6/mmm','m-3m','6/m','m-3','-1','2/m','mmm','4/m','4/mmm','-3','-3m'}; +// for i=1:numel(names) +// cs=crystalSymmetry(names{i}); +// [ad,om]=calcAngleDistribution(cs); +// oR=fundamentalRegion(cs); +// fprintf('IDX %d NAME %s MAXANGLE %.12f\n', i-1, names{i}, oR.maxAngle); +// idx=[1 25 50 75 100 125 150 175 200]; +// for k=idx +// fprintf(' P %d %.12f %.12f\n', k, om(k), ad(k)); +// end +// end +// +// The crystal-structure-index -> MTEX point-group-name map used below was +// cross-checked against EbsdLib/Core/EbsdLibConstants.h's CrystalStructure +// enum values and the ordering of LaueOps::GetAllOrientationOps(): +// 0 Hexagonal_High -> 6/mmm +// 1 Cubic_High -> m-3m +// 2 Hexagonal_Low -> 6/m +// 3 Cubic_Low -> m-3 +// 4 Triclinic -> -1 +// 5 Monoclinic -> 2/m +// 6 OrthoRhombic -> mmm +// 7 Tetragonal_Low -> 4/m +// 8 Tetragonal_High -> 4/mmm +// 9 Trigonal_Low -> -3 +// 10 Trigonal_High -> -3m + +struct SamplePoint +{ + int index; // 1-based MATLAB index into the 200-point linspace + double omega; + double ad; +}; + +void CheckDistribution(uint32_t crystalStructure, double maxAngle, const std::vector& samples) +{ + REQUIRE(random_angle_distribution::MaxMisorientationAngle(crystalStructure) == Approx(maxAngle).epsilon(1.0e-6)); + + std::vector omega(200); + for(size_t i = 0; i < omega.size(); i++) + { + omega[i] = maxAngle * static_cast(i) / static_cast(omega.size() - 1); + } + + std::vector ad = random_angle_distribution::Compute(crystalStructure, omega); + REQUIRE(ad.size() == omega.size()); + + for(const auto& sample : samples) + { + const size_t i = static_cast(sample.index - 1); + REQUIRE(omega[i] == Approx(sample.omega).epsilon(1.0e-6)); + // margin() handles the omega==maxAngle samples where the expected value + // is exactly 0.0: Approx's relative epsilon alone requires an exact + // match against 0 (scale defaults to 0), but rmag=tan(omega/2) sits + // right at its pole there, leaving ~1e-13 floating-point noise. + REQUIRE(ad[i] == Approx(sample.ad).epsilon(1.0e-6).margin(1.0e-9)); + } + + const double mean = std::accumulate(ad.cbegin(), ad.cend(), 0.0) / static_cast(ad.size()); + REQUIRE(mean == Approx(1.0).epsilon(1.0e-6)); +} +} // namespace + +TEST_CASE("RandomAngleDistribution matches MTEX for m-3m", "[RandomAngleDistribution]") +{ + CheckDistribution(ebsdlib::CrystalStructure::Cubic_High, 1.096056815241, + { + {1, 0.000000000000, 0.000000000000}, + {25, 0.132187756612, 0.073415322873}, + {50, 0.269883336416, 0.304614948103}, + {75, 0.407578916220, 0.689349878258}, + {100, 0.545274496024, 1.220337029329}, + {125, 0.682970075828, 1.887524743076}, + {150, 0.820665655632, 2.293864296709}, + {175, 0.958361235437, 1.401263125507}, + {200, 1.096056815241, 0.000000000000}, + }); +} + +TEST_CASE("RandomAngleDistribution matches MTEX for 6/mmm", "[RandomAngleDistribution]") +{ + CheckDistribution(ebsdlib::CrystalStructure::Hexagonal_High, 1.637833825000, + { + {1, 0.000000000000, 0.000000000000}, + {25, 0.197527697487, 0.122258897932}, + {50, 0.403285715703, 0.504392451654}, + {75, 0.609043733920, 0.963781985348}, + {100, 0.814801752136, 1.225761427732}, + {125, 1.020559770352, 1.436029416921}, + {150, 1.226317788568, 1.585715333366}, + {175, 1.432075806784, 1.668504346069}, + {200, 1.637833825000, 0.000000000000}, + }); +} + +TEST_CASE("RandomAngleDistribution matches MTEX for -3", "[RandomAngleDistribution]") +{ + CheckDistribution(ebsdlib::CrystalStructure::Trigonal_Low, 3.141592653590, + { + {1, 0.000000000000, 0.000000000000}, + {25, 0.378885546162, 0.213840156584}, + {50, 0.773557990080, 0.858015264909}, + {75, 1.168230433998, 1.601620985803}, + {100, 1.562902877917, 1.740726967230}, + {125, 1.957575321835, 1.612187686555}, + {150, 2.352247765753, 1.235766685169}, + {175, 2.746920209671, 0.669340528134}, + {200, 3.141592653590, 0.000000000000}, + }); +} + +TEST_CASE("RandomAngleDistribution matches MTEX for 2/m", "[RandomAngleDistribution]") +{ + CheckDistribution(ebsdlib::CrystalStructure::Monoclinic, 3.141592653590, + { + {1, 0.000000000000, 0.000000000000}, + {25, 0.378885546162, 0.142560751959}, + {50, 0.773557990080, 0.572012774922}, + {75, 1.168230433998, 1.222576203885}, + {100, 1.562902877917, 1.994223738334}, + {125, 1.957575321835, 1.861602445795}, + {150, 2.352247765753, 1.426946938453}, + {175, 2.746920209671, 0.772891378985}, + {200, 3.141592653590, 0.000000000000}, + }); +} + +TEST_CASE("RandomAngleDistribution throws on unknown crystal structure", "[RandomAngleDistribution]") +{ + const std::vector omega{0.0, 0.5, 1.0}; + REQUIRE_THROWS_AS(random_angle_distribution::Compute(ebsdlib::CrystalStructure::UnknownCrystalStructure, omega), std::invalid_argument); + REQUIRE_THROWS_AS(random_angle_distribution::MaxMisorientationAngle(ebsdlib::CrystalStructure::UnknownCrystalStructure), std::invalid_argument); +} diff --git a/Source/Test/UnitTestSupport.hpp b/Source/Test/UnitTestSupport.hpp index eda60384..b73ab58e 100644 --- a/Source/Test/UnitTestSupport.hpp +++ b/Source/Test/UnitTestSupport.hpp @@ -380,10 +380,7 @@ inline bool AlmostEqualUlpsFinal(float* A, float* B, int maxUlps) // ----------------------------------------------------------------------------- // Developer Used Macros // ----------------------------------------------------------------------------- -#define DREAM3D_TEST_FAILED(P) \ - { \ - DREAM3D_TEST_THROW_EXCEPTION(s) \ - } +#define DREAM3D_TEST_FAILED(P) {DREAM3D_TEST_THROW_EXCEPTION(s)} #define DREAM3D_REQUIRE(P) \ { \