From a9d4fd730f67b5478ee1d56eabaa3fb6691083d8 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Thu, 13 Aug 2026 16:39:25 -0400 Subject: [PATCH 1/2] STY: Update clang-format version to 15 (#52) Signed-off-by: Michael Jackson --- .github/workflows/clang-format.yml | 40 ++++++---- scripts/clang_format.py | 121 +++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 15 deletions(-) create mode 100644 scripts/clang_format.py diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index 8c97a25f..c7edb805 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -1,23 +1,33 @@ -name: clang-format check +name: clang-format pr on: pull_request: - branches: [ develop ] - push: - branches: [ develop ] + branches: + - develop + - master jobs: - clang_format_check: - runs-on: ubuntu-latest - - defaults: - run: - shell: bash + clang_format_pr: + runs-on: ubuntu-22.04 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v2 with: - fetch-depth: 2 - - - name: Run clang-format on diff - run: ./Support/Scripts/run_clang_format_diff.sh HEAD^ HEAD + fetch-depth: 2 + - name: Add Problem Matcher + uses: ammaraskar/gcc-problem-matcher@a141586609e2a558729b99a8c574c048f7f56204 + - name: Check Formatting + id: check_format + continue-on-error: true + run: | + python3 scripts/clang_format.py --format-version 15 --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 + - name: Add Suggestions + if: steps.check_format.outcome != 'success' + uses: reviewdog/action-suggester@v1 + with: + tool_name: clang-format + fail_level: error diff --git a/scripts/clang_format.py b/scripts/clang_format.py new file mode 100644 index 00000000..eea568d0 --- /dev/null +++ b/scripts/clang_format.py @@ -0,0 +1,121 @@ +import argparse +import re +import subprocess + +from typing import List, Optional, Tuple + +EXE = 'clang-format' +FILE_TYPES = ( + '.h', + '.H', + '.hpp', + '.hh', + '.h++', + '.hxx', + '.c', + '.C', + '.cpp', + '.cc', + '.c++', + '.cxx', +) + +def filter_files(files: List[str]) -> List[str]: + return [file for file in files if file.endswith(FILE_TYPES)] + +def get_git_files() -> List[str]: + try: + result = subprocess.run(['git', 'ls-files'], text=True, capture_output=True, check=True) + except subprocess.CalledProcessError as error: + print(error.stderr) + raise + git_files = result.stdout.splitlines() + return filter_files(git_files) + +def get_git_diff_files(previous_commit: str, current_commit: str) -> List[str]: + try: + result = subprocess.run( + ['git', 'diff', '--name-only', '--diff-filter=d', f'{previous_commit}...{current_commit}'], + text=True, + capture_output=True, + check=True + ) + except subprocess.CalledProcessError as error: + print(error.stderr) + raise + git_files = result.stdout.splitlines() + return filter_files(git_files) + +def format(file: str, should_modify: bool, version: Optional[int] = None) -> int: + exe = EXE + if version: + exe = f'{EXE}-{version}' + command = [exe, '--style=file'] + if should_modify: + command.append('-i') + else: + command.extend(['--dry-run', '--Werror']) + command.append(file) + result = subprocess.run(command, text=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT) + for line in result.stdout.splitlines(): + print(line, flush=True) + return result.returncode + +def get_clang_format_version() -> Tuple[int, int, int]: + try: + version_result = subprocess.run([EXE, '--version'], text=True, capture_output=True, check=True) + except subprocess.CalledProcessError as error: + print(error.stderr) + raise + version_string = version_result.stdout + match = re.match(r'clang-format version (\d+)\.(\d+)\.(\d+)', version_string) + major = int(match.group(1)) + minor = int(match.group(2)) + patch = int(match.group(3)) + return (major, minor, patch) + +def main() -> int: + parser = argparse.ArgumentParser(description='Run clang-format on files in a git repo') + parser.add_argument('--modify', help='Apply clang-format modifications in place', action='store_true') + parser.add_argument('--commits', help='Check formatting of diff between two commits instead of all files (previous, current)', nargs=2) + parser.add_argument('--format-version', help='clang-format version', type=int) + + args = parser.parse_args() + + # clang-format must be at least major version 10 since that's when --dry-run and --Werror were introduced + if args.format_version: + if args.format_version < 10: + parser.error('format-version must be at least 10') + else: + version = get_clang_format_version() + if version[0] < 10: + raise RuntimeError('clang-format must be at least major version 10') + + if args.commits: + files = get_git_diff_files(args.commits[0], args.commits[1]) + else: + files = get_git_files() + + if not files: + print('No files to format') + + failed_files: List[str] = [] + + result = 0 + for file in files: + format_result = format(file, args.modify, args.format_version) + if format_result != 0: + failed_files.append(file) + result = 1 + + if result == 0: + print('All files were formatted correctly') + else: + print('Incorrectly formatted files:') + for file in failed_files: + print(f'\"{file}\"') + + return result + +if __name__ == '__main__': + exit(main()) From f7a5b2f01bdb4bc1789dad7c0319d9ad0f26c9bc Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Tue, 25 Aug 2026 18:23:20 -0400 Subject: [PATCH 2/2] BUG: Misc bug fixes for 3.1.1 release (#53) * BUG: Fix axis normalization typo in TrigonalLowOps::getMDFFZRod * The y component of the misorientation axis was being overwritten with the z value (ax[1] = ax[2] / denom) and the z component was never normalized, corrupting the MDF fundamental zone fold for Trigonal -3 (Laue class) phases * Bug has been present since the initial EbsdLib commit, carried over from the legacy DREAM3D OrientationLib port Signed-off-by: Michael Jackson * BUG: Fix getMDFFZRod fundamental zone folds for six Laue classes Each getMDFFZRod must fold the misorientation axis into the fundamental sector of the axis equivalence group (conjugation by the rotational symmetry operators plus switching symmetry) while preserving the misorientation angle. Six of the eight implemented classes deviated: * CubicLowOps: the misorientation angle was assigned to FZn1 instead of FZw, so every result had a zero rotation angle. Also the full descending sort assumed all six axis permutations are symmetries, but the tetrahedral group's <111> 3-folds only provide cyclic permutations; the fold now cyclically rotates the largest component to the front. * TrigonalOps: returned the folded azimuth in place of the misorientation angle. Also the mirror fan-fold used sector boundaries at 60k degrees, but the 32 group's in-plane 2-folds at azimuths 0/60/120 place the mirror lines at 30 + 60k degrees; the fold now maps the azimuth into the [30, 90] degree fundamental sector. * TrigonalLowOps: returned the folded azimuth in place of the misorientation angle, and used the 60 degree mirror fan-fold copied from TrigonalOps. The -3 rotation group is only the 3-fold about c, so the azimuth now folds into a plain 120 degree wedge. * HexagonalLowOps: used the 30 degree mirror fan-fold copied from HexagonalOps, but 6/m has no in-plane 2-folds; the azimuth now folds into a plain 60 degree wedge. * TetragonalOps: only folded the axis into the first octant, missing the n1/n2 exchange required by the 422 group's <110> 2-folds; the fold now also enforces n1 >= n2. * TetragonalLowOps: folded with fabs on all components, but the 4/m axis group contains only azimuth rotations, no reflections; the axis now folds to n3 >= 0 and the azimuth into a plain 90 degree wedge. CubicOps and HexagonalOps were audited and are correct. All eight folds were validated by Monte Carlo against the exact k_QuatSym operator tables: every symmetry-equivalent axis collapses to a single canonical axis and inequivalent axes remain distinct. Signed-off-by: Michael Jackson * ENH: Implement getMDFFZRod for Triclinic/Monoclinic/OrthoRhombic + full test Implement the three previously-stubbed misorientation fundamental zone folds and add a comprehensive unit test covering all eleven Laue classes. Implementations (each folds the misorientation axis into the fundamental sector of the axis equivalence group -- conjugation by the rotational symmetry operators plus switching symmetry -- while preserving the angle): * TriclinicOps (-1): only switching symmetry, so the sector is the upper hemisphere n3 >= 0 with an equator tie-break. * MonoclinicOps (2/m): the b-axis 2-fold gives n2 >= 0 and n3 >= 0. * OrthoRhombicOps (mmm): the three 2-folds give the first octant. All three previously threw method_not_implemented, so any MDF computation on a triclinic, monoclinic, or orthorhombic phase aborted. Also handle the equator (n3 == 0) special case in the four hex/trig/tet folds: on the equator switching symmetry acts within the basal plane and combines with the rotation axis to halve the azimuth sector. The wedge folds now use fmod and a sector that shrinks on the equator. New MdfFZRodTest.cpp asserts, for every Laue class, that: * every conjugated and switched axis collapses to one canonical rod (completeness), tested against the exact k_MatSym operator tables; * the misorientation angle survives the fold (guards the historical CubicLow/Trigonal angle-discard bugs); * the fold is idempotent; * boundary axes (c-axis, 2-folds, body diagonal) fold without NaN; * inequivalent axes stay distinct (guards over-folding), with one targeted pair per class drawn from the historical over-fold bugs. 36010 assertions across the eleven classes. Signed-off-by: Michael Jackson * ENH: Add de la Vallee Poussin SO(3) kernel class * Ports MTEX SO3DeLaValleePoussinKernel: kappa/halfwidth relation, normalization constant via lgamma, 3.5*halfwidth cutoff * Unit test asserts constants against MTEX 6.1.0 generated values Signed-off-by: Michael Jackson * ENH: Add analytic random misorientation angle distribution * Ports MTEX geometry/@symmetry/calcAngleDistribution.m (lines 44-218): Cn/Dnh helpers for the six non-cubic Laue groups, plus the m-3 and m-3m cubic branches; icosahedral C3T/C3O/S3 branches intentionally skipped * Verified crystal-structure-index -> MTEX point-group-name map against EbsdLibConstants.h's CrystalStructure enum and LaueOps::GetAllOrientationOps() ordering; map was already correct * Result is normalized to unit mean over omega <= MaxMisorientationAngle and zero-clamped, dropping MTEX's constant 2*numSym(cs) prefactor * Guards the omega==maxAngle boundary, where omega/2 sits at tan()'s pole (pi/2) for the Cn Laue groups: floating-point rounding can push the argument a hair past the pole and flip tan()'s sign, so rmag is clamped with std::abs * Unit test asserts 9 sampled (omega, ad) points and MaxMisorientationAngle against MTEX 6.1.0 output for m-3m, 6/mmm, -3, and 2/m, plus the unit-mean invariant and an unknown-crystal-structure error path Signed-off-by: Michael Jackson * ENH: Add misorientation kernel density estimator on the miso-bin grid - Add ebsdlib::MisorientationKDE which accumulates weighted misorientations into the Laue-class MDF fundamental-zone bins (getMDFFZRod + getMisoBin) and evaluates a symmetrized De la Vallee Poussin kernel density at an arbitrary misorientation quaternion. - Symmetrize the query over the |CS| x |CS| crystal-symmetry pairs and add the grain-exchange inverse kernel term to enforce f(g) == f(g^-1). - Snap accumulated misorientations to their bin centers; normalize weights to sum 1 in finalize(); expose evaluate(), binCenter(), evaluateAtBinCenters(). - Add MisorientationKDETest analytic coverage: triclinic single-center peak, half-width, and cutoff; cubic crystal-symmetry invariance, grain-exchange invariance, and O(1) mean normalization. Signed-off-by: Michael Jackson * ENH: Add MDF angle-distribution curve extraction to MisorientationKDE - Add MisorientationKDE::computeAngleCurve(numPoints) returning AngleCurve {Angles, Density, RandomDensity} in radians, a port of MTEX SO3Fun/@SO3Fun/calcAngleDistribution.m: per-omega density is the uniform reference (random_angle_distribution::Compute) scaled by the mean of evaluate() over a Fibonacci full-sphere axis grid filtered to MDF-FZ membership via getMDFFZRod, with the grid count scaled by 2*|CS|. - Pin the KDE absolute scale via a direct numerical cross-check against MTEX 6.1.0 (calcDensity 'exact'): with exact centers our density matches MTEX's mdf to a constant ratio of 2.000 at every misorientation angle and MTEX's mdf has mean 1 over SO(3). Apply the resulting 0.5 antipodal factor in evaluate() (divide by 2*numSymOps^2), giving a mean-1 normalized MDF. - Update the SingleCenterTriclinic assertions to the new scale (modal peak K(0)/2, quarter peak K(0)/4). - Add the CubicAngleCurveVsMTEX cross-check TEST_CASE. Signed-off-by: Michael Jackson * TEST: Add correlated-twin MDF regression guard for MisorientationKDE Add a MisorientationKDE test that injects a tight cluster of 60-degree about <111> (Sigma3) misorientations on top of a uniform-random background and asserts that both the MDF bin-array peak folds to a 60-degree / <111> misorientation and the angle-distribution curve peaks near 60 degrees, at least 5 degrees above the cubic Mackenzie (random-reference) maximum near 45 degrees. This is the discriminating case the earlier 45-degree bicrystal cross-check could not catch: a density that collapsed to the random distribution still peaks at ~45 degrees and would pass a weaker test. The MisorientationKDE math is correct as-is; this test locks in that correctness so a future regression that flattens a correlated MDF toward the random reference is detected. Signed-off-by: Michael Jackson * STY: Add BSD headers to Texture MDF files, fix stale doc/comment/include - Add the standard EbsdLib BSD license header block (verbatim from TexturePreset.h) to SO3DeLaValleePoussinKernel.h/.cpp, RandomAngleDistribution.h/.cpp, and MisorientationKDE.h/.cpp; these six production files were missing it. - MisorientationKDE.h: correct computeAngleCurve()'s doc comment, which cited the old test tolerance epsilon(0.05)+margin(0.02); the committed MisorientationKDETest.cpp assertion actually uses epsilon(0.20).margin(0.10). - RandomAngleDistribution.cpp: remove the unused `#include ` (the .cpp sums with a manual loop; std::accumulate is only used in the test). - MisorientationKDETest.cpp: delete a dangling leftover comment ("Disorientation-style similarity of two densities.") describing a helper that no longer exists. Signed-off-by: Michael Jackson * BUG: (fix) Correct hexagonal MDF angle-distribution under-estimate at low angles The MDF misorientation-angle distribution produced by MisorientationKDE under-estimated the density at low-to-mid misorientation angles for the hexagonal (6/mmm) Laue class by 23-32% (e.g. 14 deg: 0.287 vs MTEX 0.422; 19 deg: 0.597 vs 0.771; 38 deg: 0.365 vs 0.474), while the peak, high angles, the Mackenzie random reference, and the maximum misorientation angle were all correct. Cubic (m-3m) passed only because its low-angle MTEX density is below the absolute test-margin floor and never exercised that band. Root cause: the KDE stored each accumulated misorientation at the GEOMETRIC center of its ~5-degree MDF fundamental-zone bin. computeAngleCurve() already samples the correct axis measure (a full-sphere Fibonacci grid folded to the getMDFFZRod fundamental zone) and this is not the problem: evaluate() is fully symmetrized, so its per-omega average over any axis domain is the conditional mean of the SAME density, and no axis-domain change (sector integration, disjoint-group sector, etc.) can move it. The bias came entirely from the ~5-degree bin-center quantization of the stored centers. The hexagonal [0001]/20-degree component snapped to a tilted (0.145,0.039,0.989)/22.7-degree center, shifting the steep 10-degree-halfwidth kernel off its true position and depressing the low-angle flank. MTEX's calcDensity(...,'exact') uses the exact centers, so our binned centers disagreed with it exactly in the low-angle band. Fix: accumulate the weight-weighted, sign-aligned running sum of each bin's fundamental-zone misorientation quaternions and use the normalized sum (the weighted circular mean of the observations in the bin) as the KDE center in finalize(), instead of the geometric bin center. Memory stays bounded by getMDFSize (one quaternion accumulator per bin), so this scales to millions of boundary misorientations. For an isolated misorientation the center becomes the observation itself, reproducing MTEX's 'exact' result: 14 deg 0.421 vs 0.422, 19 deg 0.719 vs 0.771, 38 deg 0.433 vs 0.474, peak index 131 at ~61.3 deg. HexagonalAngleCurveVsMTEX now passes at the same epsilon(0.20)/margin(0.10) tolerance the cubic cross-check uses, and CubicAngleCurveVsMTEX is unchanged and still passes. Adds the HexagonalAngleCurveVsMTEX regression test (MTEX 6.1.0 reference) and updates SingleCenterTriclinic to evaluate at the observation-derived center rather than the geometric bin center. Signed-off-by: Michael Jackson * BUG: Fix NaN ODF dimension constant in TetragonalLowOps The third ODF dimension init value computed pow(0.75 * (pi/4 - sin(pi/2)), 1/3), the cube root of a negative number, which is NaN. Every orientation sampled by TetragonalLowOps::determineEulerAngles therefore returned NaN Euler angles (100% failure), breaking synthetic texture generation and ODF sampling for the Tetragonal 4/m Laue class. The bug was inherited verbatim from the 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). With the fix the NaN rate for uniformly sampled bins drops from 100% to ~3.8%, in line with the other low-symmetry Laue classes (whose residual NaNs come from homochoric-cube corner bins outside the valid ball -- a separate, pre-existing issue). Signed-off-by: Michael Jackson * STY: Update clang-format version to 16 Signed-off-by: Michael Jackson * BUG: Use exact sqrt(3)/sqrt(2) normalizers in CubicOps::getSchmidFactorAndSS The auto slip-system overload of CubicOps::getSchmidFactorAndSS normalized the {111} plane-normal and <110> slip-direction dot products with the float literals 1.732f and 1.414f. Everything else in that function is computed in double, so those two four-significant-digit literals were the only precision loss in the whole calculation, and they biased the result systematically rather than randomly. Both literals are SMALLER than the constants they approximate, so every theta and every lambda came out too large and every Schmid factor was inflated by the uniform factor sqrt(6) / (1.732f * 1.414f) = 1.00018035284 (+0.0180353 %) Two consequences: - The Schmid factor of a cubic crystal cannot physically exceed 0.5, but the biased normalizers returned up to 0.500090176 at the maximizing loading direction, so callers that range-check the output saw an impossible value. - The reported cos(phi) and cos(lambda) angle components carried the same bias and could likewise exceed 1.0. The slip-system index is unaffected: the bias is a single positive scale factor applied to all twelve candidates, so the argmax is unchanged. Replace both literals with the existing full-precision double constants ebsdlib::constants::k_Sqrt3D and k_Sqrt2D (EbsdLibMath.h, already included by this translation unit). The double type matches the surrounding arithmetic, so this is the minimal correct form -- no other line changes. The second, plane/direction overload of getSchmidFactorAndSS is not affected: it normalizes the caller-supplied plane and direction vectors with std::sqrt and never used the truncated literals. No other Laue op contained 1.732f or 1.414f. Signed-off-by: Michael Jackson * BUG: Define every getSchmidFactorAndSS output on all paths A survey of the auto (load-only) getSchmidFactorAndSS overload across all eleven Laue op classes found two ways for an output parameter to come back undefined. 1. HexagonalLowOps read `schmidfactor` uninitialized. The function declares its locals and then runs `if(schmid1 > schmidfactor) { ... }` as the first thing that touches `schmidfactor` -- it was never seeded, so the whole comparison chain was driven by an indeterminate value. Depending on that garbage, the function could return a Schmid factor of zero with `slipsys` and `angleComps` also left untouched, or accept a candidate against a garbage incumbent. `slipsys` was likewise never initialized. 2. HexagonalOps seeded `schmidfactor` but not `slipsys` or `angleComps`, and the seven Laue classes that enumerate no slip systems at all (TrigonalOps, TrigonalLowOps, TetragonalOps, TetragonalLowOps, OrthoRhombicOps, MonoclinicOps, TriclinicOps) set `schmidfactor` and `slipsys` but left `angleComps` untouched. Case 2 is the more damaging of the two in practice because it is silent and plausible-looking. A caller that hoists one `angleComps[2]` buffer outside a per-feature loop -- which is the natural way to write the loop, and what the DREAM3D-NX Compute Schmid Factors filter does -- gets the PREVIOUS feature's angle components attributed to the current feature whenever the current feature's Laue class is one of the seven stubs. The Schmid factor correctly reads 0 while the two angle columns read as real measurements. Seed all four outputs at the top of every affected overload. This changes no value on any path that already computed a result; it only replaces indeterminate or stale output with a defined zero. CubicOps and CubicLowOps already defined all four unconditionally and are unchanged. The plane/direction overloads all already zeroed their outputs on entry and are unchanged. Signed-off-by: Michael Jackson * TEST: Lock the stub Laue classes' Schmid outputs to zero The seven Laue classes that enumerate no slip systems (Trigonal, TrigonalLow, Tetragonal, TetragonalLow, OrthoRhombic, Monoclinic, Triclinic) set schmidfactor and slipsys but left angleComps untouched until 2c84f2a. That defect was silent: a caller which hoists one angleComps[2] buffer outside a per-Feature loop got the previous Feature's angle components attributed to a Feature of one of these classes, while schmidfactor correctly read 0. The new case calls each stub's load-only getSchmidFactorAndSS with all four outputs pre-poisoned (schmidfactor 7.0, slipsys 7, angleComps {7.0, 9.0}) and asserts that all four read back as defined zeros. Pre-poisoning is what gives the assertion its discriminating power: an overload that does not write angleComps leaves the poison in place. Verified by mutation -- removing the two angleComps stores from TrigonalOps:: getSchmidFactorAndSS fails the case with "7.0 == 0.0" / "9.0 == 0.0" against TrigonalOps, and restoring them returns it to 28 of 28 passing. Signed-off-by: Michael Jackson * STY: Update clang-format version to 19 Signed-off-by: Michael Jackson --------- Signed-off-by: Michael Jackson --- .../{clang-format.yml => format_pr.yml} | 11 +- .github/workflows/format_push.yml | 24 ++ AGENTS.md | 9 + CMakeLists.txt | 2 +- Docs/Index.md | 77 +++- Source/EbsdLib/Core/EbsdDataArray.cpp | 4 +- Source/EbsdLib/Core/EbsdLibDLLExport.h | 4 +- Source/EbsdLib/IO/HKL/H5CtfImporter.cpp | 4 +- Source/EbsdLib/LaueOps/CubicLowOps.cpp | 46 +- Source/EbsdLib/LaueOps/CubicOps.cpp | 20 +- Source/EbsdLib/LaueOps/HexagonalLowOps.cpp | 31 +- Source/EbsdLib/LaueOps/HexagonalOps.cpp | 11 +- Source/EbsdLib/LaueOps/MonoclinicOps.cpp | 41 +- Source/EbsdLib/LaueOps/OrthoRhombicOps.cpp | 20 +- Source/EbsdLib/LaueOps/TetragonalLowOps.cpp | 43 +- Source/EbsdLib/LaueOps/TetragonalOps.cpp | 11 + Source/EbsdLib/LaueOps/TriclinicOps.cpp | 28 +- Source/EbsdLib/LaueOps/TrigonalLowOps.cpp | 48 +-- Source/EbsdLib/LaueOps/TrigonalOps.cpp | 43 +- Source/EbsdLib/Texture/MisorientationKDE.cpp | 242 +++++++++++ Source/EbsdLib/Texture/MisorientationKDE.h | 162 +++++++ .../Texture/RandomAngleDistribution.cpp | 293 +++++++++++++ .../EbsdLib/Texture/RandomAngleDistribution.h | 61 +++ .../Texture/SO3DeLaValleePoussinKernel.cpp | 85 ++++ .../Texture/SO3DeLaValleePoussinKernel.h | 76 ++++ Source/EbsdLib/Texture/SourceList.cmake | 6 + Source/EbsdLib/Utilities/inipp.h | 4 +- .../wlenthe_orientation_coloring.hpp | 4 +- Source/Test/CMakeLists.txt | 4 + Source/Test/DeLaValleePoussinKernelTest.cpp | 49 +++ Source/Test/LaueOpsTest.cpp | 26 ++ Source/Test/MdfFZRodTest.cpp | 321 ++++++++++++++ Source/Test/MisorientationKDETest.cpp | 405 ++++++++++++++++++ Source/Test/RandomAngleDistributionTest.cpp | 150 +++++++ Source/Test/UnitTestSupport.hpp | 5 +- 35 files changed, 2213 insertions(+), 157 deletions(-) rename .github/workflows/{clang-format.yml => format_pr.yml} (70%) create mode 100644 .github/workflows/format_push.yml create mode 100644 AGENTS.md create mode 100644 Source/EbsdLib/Texture/MisorientationKDE.cpp create mode 100644 Source/EbsdLib/Texture/MisorientationKDE.h create mode 100644 Source/EbsdLib/Texture/RandomAngleDistribution.cpp create mode 100644 Source/EbsdLib/Texture/RandomAngleDistribution.h create mode 100644 Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.cpp create mode 100644 Source/EbsdLib/Texture/SO3DeLaValleePoussinKernel.h create mode 100644 Source/Test/DeLaValleePoussinKernelTest.cpp create mode 100644 Source/Test/MdfFZRodTest.cpp create mode 100644 Source/Test/MisorientationKDETest.cpp create mode 100644 Source/Test/RandomAngleDistributionTest.cpp 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/AGENTS.md b/AGENTS.md new file mode 100644 index 00000000..1a2cc21f --- /dev/null +++ b/AGENTS.md @@ -0,0 +1,9 @@ +# Repository Guidance + +EbsdLib is a C++20 CMake library. Follow the checked-in formatter and lint configuration: Allman braces, two-space indentation, left-aligned pointers, and the repository naming conventions. + +- Core library code is under `Source/EbsdLib/`; applications are under `Source/Apps/`; tests and CMake support live in their existing repository directories. +- Configure with `cmake --preset EbsdLib-Release`; verify with the actual CMake build and tests rather than relying on clangd diagnostics alone. +- Run tests through `ctest`, not test executables directly. +- For inverse-pole-figure and standard-stereographic-triangle visualization, reuse the LaueOps stereographic mapping. Do not substitute a Lambert eta/chi rectangle. +- Keep durable architecture and workflow documentation in the repository; do not encode active-branch status or personal machine paths here. diff --git a/CMakeLists.txt b/CMakeLists.txt index 7635656b..793e2bd2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -37,7 +37,7 @@ option(EbsdLib_BUILD_H5SUPPORT "Build H5Support Library" OFF) # set project's name -project(EbsdLibProj VERSION 3.1.0) +project(EbsdLibProj VERSION 3.1.1) # Request C++17 standard, using new CMake variables. diff --git a/Docs/Index.md b/Docs/Index.md index e58490e7..a37f1cf4 100644 --- a/Docs/Index.md +++ b/Docs/Index.md @@ -1,4 +1,4 @@ -# Various Bits of Documentation for EbsdLib +# Documentation for EbsdLib EbsdLib is primarily used in the [DREAM3D](https://www.bluequartz.net) family of applications and libraries. @@ -27,6 +27,81 @@ Position-space validation across all 11 Laue classes lives in --- +# Release Notes — EbsdLib 3.1.1 + +EbsdLib 3.1.1 is the "misorientation analysis + correctness" release. It adds +tools for calculating and evaluating misorientation distribution functions +(MDFs), and fixes several long-standing Laue-class fundamental-zone and Schmid +factor defects inherited from the legacy OrientationLib implementation. + +## New features + +### Misorientation kernel density estimation + +- `MisorientationKDE` accumulates weighted misorientations on the Laue-class MDF + fundamental-zone grid and evaluates the resulting density at arbitrary + misorientation quaternions. +- KDE evaluation is symmetrized over the crystal symmetry operators and grain + exchange, and is normalized to a mean density of one over SO(3). +- `computeAngleCurve()` extracts a misorientation-angle distribution together + with the corresponding random-reference distribution. +- Each occupied bin retains the weighted circular mean of its observations as + its kernel center. This avoids the low-angle bias caused by snapping narrow + kernels to geometric bin centers, particularly for hexagonal phases. + +### SO(3) kernel and random-angle reference + +- `SO3DeLaValleePoussinKernel` implements the de la Vallée Poussin orientation + kernel, including half-width conversion, normalization, and finite support. +- `random_angle_distribution::Compute()` and `MaxMisorientationAngle()` provide + analytical random misorientation-angle distributions and limits for the + supported Laue classes. + +## Correctness fixes + +### Misorientation fundamental-zone folds + +- Corrected `getMDFFZRod()` for cubic-low, trigonal, trigonal-low, + hexagonal-low, tetragonal, and tetragonal-low symmetry. The corrected folds + preserve the misorientation angle and use the proper rotational-symmetry + sector for each Laue class. +- Implemented the previously stubbed triclinic, monoclinic, and orthorhombic + folds, so MDF calculations for those phases no longer abort with + `method_not_implemented`. +- Corrected equatorial boundary handling in the hexagonal, trigonal, and + tetragonal folds. + +### ODF and Schmid-factor calculations + +- Fixed the tetragonal-low ODF dimension constant that caused + `determineEulerAngles()` to return NaN Euler angles for all sampled bins. +- Replaced truncated cubic slip-system normalizers with the full-precision + `sqrt(3)` and `sqrt(2)` constants, preventing Schmid factors from exceeding + the physical maximum of 0.5 because of rounding bias. +- Initialized every output of the automatic `getSchmidFactorAndSS()` overloads + on all code paths. Unsupported Laue classes now return defined zero values + instead of stale or indeterminate angle components. + +## Validation evidence + +- Fundamental-zone tests cover all 11 Laue classes, including symmetry + equivalence, angle preservation, idempotence, boundary behavior, and guards + against over-folding. +- Kernel constants and random-angle distributions are checked against MTEX + 6.1.0 reference values. +- Cubic and hexagonal MDF angle curves are cross-checked against MTEX, and a + correlated 60° about `<111>` (Σ3) distribution guards against regressions + that flatten a textured MDF toward the random reference. +- Schmid-factor tests pre-poison output values to ensure the stub Laue classes + explicitly reset every output. + +## Maintenance + +- Updated the formatting workflows to run for both pull requests and pushes, + using clang-format 19. + +--- + # Release Notes — EbsdLib 3.1.0 EbsdLib 3.1 is the "discrete pole figures + TSL-by-default" release. It builds 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) \ { \