BUG: Misc bug fixes for 3.1.1 release (#53)

[imikejackson]

• 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

• 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.

• 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.

• 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

• 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

• 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.

• 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.

• 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.

• 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 <numeric> (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.

• 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.

• 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).

• STY: Update clang-format version to 16

• 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.

• 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.

• 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.

• STY: Update clang-format version to 19

---