From b0d864f474b14f758bd8cede182f4e5ead8ca1c6 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 22:50:06 -0400 Subject: [PATCH 01/16] feat: add ColorSpaceUtils with HSL/HSV/RGB conversions Header-only utility in ebsdlib::color namespace providing hslToRgb, hslToHsv, and rgbToBytes functions, with Catch2 unit tests covering primary colors, achromatic cases, and round-trip value bounds. Co-Authored-By: Claude Sonnet 4.6 --- Source/EbsdLib/Utilities/ColorSpaceUtils.hpp | 94 ++++++++++++++++++++ Source/EbsdLib/Utilities/SourceList.cmake | 1 + Source/Test/CMakeLists.txt | 1 + Source/Test/ColorSpaceUtilsTest.cpp | 89 ++++++++++++++++++ 4 files changed, 185 insertions(+) create mode 100644 Source/EbsdLib/Utilities/ColorSpaceUtils.hpp create mode 100644 Source/Test/ColorSpaceUtilsTest.cpp diff --git a/Source/EbsdLib/Utilities/ColorSpaceUtils.hpp b/Source/EbsdLib/Utilities/ColorSpaceUtils.hpp new file mode 100644 index 0000000..b70c6bd --- /dev/null +++ b/Source/EbsdLib/Utilities/ColorSpaceUtils.hpp @@ -0,0 +1,94 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" + +#include +#include +#include +#include + +namespace ebsdlib +{ +namespace color +{ + +/** + * @brief Convert HSL to RGB. All inputs and outputs in [0, 1]. + * @param h Hue in [0, 1) where 0=red, 1/3=green, 2/3=blue + * @param s Saturation in [0, 1] + * @param l Lightness in [0, 1] + * @return {r, g, b} each in [0, 1] + */ +inline std::array hslToRgb(double h, double s, double l) +{ + double c = (1.0 - std::abs(2.0 * l - 1.0)) * s; + double hp = h * 6.0; + double x = c * (1.0 - std::abs(std::fmod(hp, 2.0) - 1.0)); + double m = l - c / 2.0; + + double r1 = 0.0; + double g1 = 0.0; + double b1 = 0.0; + + if(hp < 1.0) + { + r1 = c; + g1 = x; + b1 = 0.0; + } + else if(hp < 2.0) + { + r1 = x; + g1 = c; + b1 = 0.0; + } + else if(hp < 3.0) + { + r1 = 0.0; + g1 = c; + b1 = x; + } + else if(hp < 4.0) + { + r1 = 0.0; + g1 = x; + b1 = c; + } + else if(hp < 5.0) + { + r1 = x; + g1 = 0.0; + b1 = c; + } + else + { + r1 = c; + g1 = 0.0; + b1 = x; + } + + return {std::clamp(r1 + m, 0.0, 1.0), std::clamp(g1 + m, 0.0, 1.0), std::clamp(b1 + m, 0.0, 1.0)}; +} + +/** + * @brief Convert HSL to HSV. All inputs and outputs in [0, 1]. + */ +inline std::array hslToHsv(double h, double s, double l) +{ + double l2 = 2.0 * l; + double s2 = s * ((l2 <= 1.0) ? l2 : (2.0 - l2)); + double v = (l2 + s2) / 2.0; + double sv = (l2 + s2 > 1e-12) ? (2.0 * s2 / (l2 + s2)) : 0.0; + return {h, sv, v}; +} + +/** + * @brief Convert RGB [0,1] to 8-bit [0,255] clamped. + */ +inline std::array rgbToBytes(double r, double g, double b) +{ + return {static_cast(std::clamp(r * 255.0, 0.0, 255.0)), static_cast(std::clamp(g * 255.0, 0.0, 255.0)), static_cast(std::clamp(b * 255.0, 0.0, 255.0))}; +} + +} // namespace color +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index 0e7fc72..fe772f5 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -23,6 +23,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/CanvasUtilities.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureCompositor.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/inipp.h + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ColorSpaceUtils.hpp ) set(EbsdLib_${DIR_NAME}_SRCS diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index c0754d1..c109d24 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -54,6 +54,7 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/DirectionalStatsTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/UnitTestCommon.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/UnitTestCommon.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/ColorSpaceUtilsTest.cpp ) diff --git a/Source/Test/ColorSpaceUtilsTest.cpp b/Source/Test/ColorSpaceUtilsTest.cpp new file mode 100644 index 0000000..c6a5839 --- /dev/null +++ b/Source/Test/ColorSpaceUtilsTest.cpp @@ -0,0 +1,89 @@ +#include + +#include "EbsdLib/Utilities/ColorSpaceUtils.hpp" + +TEST_CASE("ebsdlib::ColorSpaceUtils::HslToRgb", "[EbsdLib][ColorSpaceUtils]") +{ + SECTION("Pure Red") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(0.0, 1.0, 0.5); + REQUIRE(r == Approx(1.0).margin(1e-6)); + REQUIRE(g == Approx(0.0).margin(1e-6)); + REQUIRE(b == Approx(0.0).margin(1e-6)); + } + SECTION("Pure Green") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(1.0 / 3.0, 1.0, 0.5); + REQUIRE(r == Approx(0.0).margin(1e-6)); + REQUIRE(g == Approx(1.0).margin(1e-6)); + REQUIRE(b == Approx(0.0).margin(1e-6)); + } + SECTION("Pure Blue") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(2.0 / 3.0, 1.0, 0.5); + REQUIRE(r == Approx(0.0).margin(1e-6)); + REQUIRE(g == Approx(0.0).margin(1e-6)); + REQUIRE(b == Approx(1.0).margin(1e-6)); + } + SECTION("White") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(0.0, 0.0, 1.0); + REQUIRE(r == Approx(1.0).margin(1e-6)); + REQUIRE(g == Approx(1.0).margin(1e-6)); + REQUIRE(b == Approx(1.0).margin(1e-6)); + } + SECTION("Black") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(0.0, 0.0, 0.0); + REQUIRE(r == Approx(0.0).margin(1e-6)); + REQUIRE(g == Approx(0.0).margin(1e-6)); + REQUIRE(b == Approx(0.0).margin(1e-6)); + } + SECTION("50% Gray") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(0.0, 0.0, 0.5); + REQUIRE(r == Approx(0.5).margin(1e-6)); + REQUIRE(g == Approx(0.5).margin(1e-6)); + REQUIRE(b == Approx(0.5).margin(1e-6)); + } + SECTION("Yellow (H=60deg)") + { + auto [r, g, b] = ebsdlib::color::hslToRgb(1.0 / 6.0, 1.0, 0.5); + REQUIRE(r == Approx(1.0).margin(1e-6)); + REQUIRE(g == Approx(1.0).margin(1e-6)); + REQUIRE(b == Approx(0.0).margin(1e-6)); + } +} + +TEST_CASE("ebsdlib::ColorSpaceUtils::HslToHsv", "[EbsdLib][ColorSpaceUtils]") +{ + SECTION("Full saturation, mid lightness -> V=1, S=1") + { + auto [h, s, v] = ebsdlib::color::hslToHsv(0.0, 1.0, 0.5); + REQUIRE(h == Approx(0.0)); + REQUIRE(s == Approx(1.0)); + REQUIRE(v == Approx(1.0)); + } + SECTION("Zero saturation -> S_hsv = 0") + { + auto [h, s, v] = ebsdlib::color::hslToHsv(0.5, 0.0, 0.5); + REQUIRE(s == Approx(0.0)); + REQUIRE(v == Approx(0.5)); + } +} + +TEST_CASE("ebsdlib::ColorSpaceUtils::RoundTrip", "[EbsdLib][ColorSpaceUtils]") +{ + for(double hue = 0.0; hue < 1.0; hue += 0.1) + { + auto [r, g, b] = ebsdlib::color::hslToRgb(hue, 1.0, 0.5); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + double maxVal = std::max({r, g, b}); + REQUIRE(maxVal == Approx(1.0).margin(1e-6)); + } +} From e59e42a2d1c83ded737501199280c8652ab60e5a Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 22:51:20 -0400 Subject: [PATCH 02/16] feat: add IColorKey abstract interface for pluggable IPF color strategies Co-Authored-By: Claude Sonnet 4.6 --- Source/EbsdLib/Utilities/IColorKey.hpp | 63 +++++++++++++++++++++++ Source/EbsdLib/Utilities/SourceList.cmake | 1 + 2 files changed, 64 insertions(+) create mode 100644 Source/EbsdLib/Utilities/IColorKey.hpp diff --git a/Source/EbsdLib/Utilities/IColorKey.hpp b/Source/EbsdLib/Utilities/IColorKey.hpp new file mode 100644 index 0000000..c19063e --- /dev/null +++ b/Source/EbsdLib/Utilities/IColorKey.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" + +#include +#include +#include +#include + +namespace ebsdlib +{ + +/** + * @brief Abstract interface for IPF color key strategies. + * + * Maps a crystal direction (already projected into the fundamental sector) + * to an RGB color. Implementations include TSL (traditional) and + * Nolze-Hielscher (perceptually improved). + * + * Two overloads are provided: + * - direction2Color(Vec3): takes a 3D unit direction vector (preferred for N-H) + * - direction2Color(eta, chi, angleLimits): takes spherical coords (TSL compatibility) + */ +class EbsdLib_EXPORT IColorKey +{ +public: + using Pointer = std::shared_ptr; + using Vec3 = std::array; + + virtual ~IColorKey() = default; + + /** + * @brief Map a unit crystal direction vector (in the fundamental sector) to an RGB color. + * This is the primary interface. The direction must already be projected into the SST. + * @param direction Unit direction vector {x, y, z} in the fundamental sector + * @return {R, G, B} each in [0.0, 1.0] + */ + virtual Vec3 direction2Color(const Vec3& direction) const = 0; + + /** + * @brief Map a crystal direction via spherical coordinates to an RGB color. + * Provided for backward compatibility with the TSL pipeline. + * Default implementation converts to a direction vector and calls the Vec3 overload. + * @param eta Azimuthal angle of the direction (radians) + * @param chi Polar angle of the direction from z-axis (radians) + * @param angleLimits {etaMin, etaMax, chiMax} from the LaueOps subclass + * @return {R, G, B} each in [0.0, 1.0] + */ + virtual Vec3 direction2Color(double eta, double chi, const Vec3& angleLimits) const + { + // Default: convert spherical to Cartesian and delegate + double sinChi = std::sin(chi); + Vec3 dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + return direction2Color(dir); + } + + /** + * @brief Human-readable name of this color key. + */ + virtual std::string name() const = 0; +}; + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index fe772f5..2990909 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -24,6 +24,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureCompositor.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/inipp.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ColorSpaceUtils.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/IColorKey.hpp ) set(EbsdLib_${DIR_NAME}_SRCS From bf6803126bb2053d0082d69896e7b989ef07b83b Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:05:17 -0400 Subject: [PATCH 03/16] feat: add FundamentalSectorGeometry with polar coordinate computation for all 11 Laue groups Implements sector geometry definitions (boundary normals, vertices, barycenter) and normalized polar coordinate computation (radius, azimuth) for all 11 Laue groups. Each sector is defined by half-space boundary normals with factory methods for cubicHigh (m-3m), cubicLow (m-3), hexagonalHigh (6/mmm), hexagonalLow (6/m), tetragonalHigh (4/mmm), tetragonalLow (4/m), trigonalHigh (-3m), trigonalLow (-3), orthorhombic (mmm), monoclinic (2/m), and triclinic (-1). Includes 8 Catch2 unit tests covering vertex geometry, polar coordinates, isInside checks, barycenter validity, and azimuthal correction. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/FundamentalSectorGeometry.cpp | 657 ++++++++++++++++++ .../Utilities/FundamentalSectorGeometry.hpp | 66 ++ Source/EbsdLib/Utilities/SourceList.cmake | 2 + Source/Test/CMakeLists.txt | 1 + Source/Test/FundamentalSectorGeometryTest.cpp | 286 ++++++++ 5 files changed, 1012 insertions(+) create mode 100644 Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp create mode 100644 Source/EbsdLib/Utilities/FundamentalSectorGeometry.hpp create mode 100644 Source/Test/FundamentalSectorGeometryTest.cpp diff --git a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp new file mode 100644 index 0000000..6924133 --- /dev/null +++ b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp @@ -0,0 +1,657 @@ +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" + +#include +#include +#include +#include +#include + +namespace ebsdlib +{ + +// ----------------------------------------------------------------------- +// Constructor +// ----------------------------------------------------------------------- +FundamentalSectorGeometry::FundamentalSectorGeometry(std::vector boundaryNormals, std::vector vertices, + std::string colorKeyMode, int32_t supergroupIndex) +: m_BoundaryNormals(std::move(boundaryNormals)) +, m_Vertices(std::move(vertices)) +, m_ColorKeyMode(std::move(colorKeyMode)) +, m_SupergroupIndex(supergroupIndex) +{ + computeBarycenter(); + precomputeAzimuthalCorrection(); +} + +// ----------------------------------------------------------------------- +// Accessors +// ----------------------------------------------------------------------- +const FundamentalSectorGeometry::Vec3& FundamentalSectorGeometry::barycenter() const +{ + return m_Barycenter; +} + +const std::vector& FundamentalSectorGeometry::vertices() const +{ + return m_Vertices; +} + +const std::vector& FundamentalSectorGeometry::boundaryNormals() const +{ + return m_BoundaryNormals; +} + +const std::string& FundamentalSectorGeometry::colorKeyMode() const +{ + return m_ColorKeyMode; +} + +int32_t FundamentalSectorGeometry::supergroupIndex() const +{ + return m_SupergroupIndex; +} + +// ----------------------------------------------------------------------- +// Vector math helpers +// ----------------------------------------------------------------------- +FundamentalSectorGeometry::Vec3 FundamentalSectorGeometry::vecNormalize(const Vec3& v) +{ + double len = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + if(len < 1.0e-15) + { + return {0.0, 0.0, 0.0}; + } + return {v[0] / len, v[1] / len, v[2] / len}; +} + +FundamentalSectorGeometry::Vec3 FundamentalSectorGeometry::vecCross(const Vec3& a, const Vec3& b) +{ + return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]}; +} + +double FundamentalSectorGeometry::vecDot(const Vec3& a, const Vec3& b) +{ + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +double FundamentalSectorGeometry::vecAngle(const Vec3& a, const Vec3& b) +{ + double d = vecDot(a, b); + d = std::clamp(d, -1.0, 1.0); + return std::acos(d); +} + +FundamentalSectorGeometry::Vec3 FundamentalSectorGeometry::vecNeg(const Vec3& v) +{ + return {-v[0], -v[1], -v[2]}; +} + +// ----------------------------------------------------------------------- +// computeBarycenter +// ----------------------------------------------------------------------- +void FundamentalSectorGeometry::computeBarycenter() +{ + if(m_Vertices.empty()) + { + // For triclinic: use the north pole as the center of the hemisphere + m_Barycenter = {0.0, 0.0, 1.0}; + return; + } + + Vec3 sum = {0.0, 0.0, 0.0}; + for(const auto& v : m_Vertices) + { + sum[0] += v[0]; + sum[1] += v[1]; + sum[2] += v[2]; + } + m_Barycenter = vecNormalize(sum); +} + +// ----------------------------------------------------------------------- +// isInside +// ----------------------------------------------------------------------- +bool FundamentalSectorGeometry::isInside(const Vec3& h) const +{ + for(const auto& normal : m_BoundaryNormals) + { + if(vecDot(h, normal) < -1.0e-10) + { + return false; + } + } + return true; +} + +// ----------------------------------------------------------------------- +// polarCoordinates +// ----------------------------------------------------------------------- +std::pair FundamentalSectorGeometry::polarCoordinates(const Vec3& h) const +{ + constexpr double k_Pi = 3.14159265358979323846; + + // Singularity guard: if h is at the barycenter + double angleToCenter = vecAngle(h, m_Barycenter); + if(angleToCenter < 1.0e-10) + { + return {0.0, 0.0}; + } + + // ------------------------------------------------------------------- + // RADIUS: normalized distance from center to boundary + // ------------------------------------------------------------------- + // Normal to the great circle through center and h + Vec3 gcNormal = vecNormalize(vecCross(m_Barycenter, h)); + + double radius = std::numeric_limits::infinity(); + double distCenterH = angleToCenter; // angle from center to h + + for(const auto& normal : m_BoundaryNormals) + { + // Intersection of the great circle (center,h) with boundary plane N_j. + // The two great circles intersect at: P = normalize(cross(N_j, gcNormal)) and -P. + Vec3 bp = vecNormalize(vecCross(normal, gcNormal)); + + // Choose the intersection point on the same side as h relative to center. + // We want the point P such that the arc center->P passes through or near h. + // Test: if dot(h, bp) < 0, flip to -bp (choose the hemisphere containing h). + if(vecDot(h, bp) < 0.0) + { + bp = vecNeg(bp); + } + + double distCenterBp = vecAngle(m_Barycenter, bp); + + double ratio; + if(distCenterBp < 1.0e-10) + { + // Center is on this boundary -- boundary is at distance 0, but the sector + // is on the positive side. This shouldn't happen for well-defined sectors. + // Treat as no constraint. + continue; + } + else + { + ratio = distCenterH / distCenterBp; + } + + // Handle NaN + if(std::isnan(ratio)) + { + ratio = 1.0; + } + + radius = std::min(radius, ratio); + } + + if(std::isinf(radius)) + { + radius = 0.0; + } + radius = std::clamp(radius, 0.0, 1.0); + + // ------------------------------------------------------------------- + // AZIMUTHAL ANGLE: angle in the tangent plane at the barycenter + // ------------------------------------------------------------------- + // Reference direction: project z-axis onto tangent plane at barycenter + Vec3 ref = {0.0, 0.0, 1.0}; + // If barycenter is very close to z-axis, use x-axis instead + if(std::abs(vecDot(ref, m_Barycenter)) > 0.99) + { + ref = {1.0, 0.0, 0.0}; + } + + // Project ref onto tangent plane at barycenter: rx = ref - dot(ref, center) * center + double refDotCenter = vecDot(ref, m_Barycenter); + Vec3 rx = {ref[0] - refDotCenter * m_Barycenter[0], ref[1] - refDotCenter * m_Barycenter[1], ref[2] - refDotCenter * m_Barycenter[2]}; + rx = vecNormalize(rx); + + // ry = center x rx (right-hand rule in tangent plane) + Vec3 ry = vecNormalize(vecCross(m_Barycenter, rx)); + + // Direction from center to h in tangent plane (not normalized -- fine for atan2) + double hDotCenter = vecDot(h, m_Barycenter); + Vec3 dv = {h[0] - hDotCenter * m_Barycenter[0], h[1] - hDotCenter * m_Barycenter[1], h[2] - hDotCenter * m_Barycenter[2]}; + + double rho = std::atan2(vecDot(ry, dv), vecDot(rx, dv)); + rho = std::fmod(rho + 2.0 * k_Pi, 2.0 * k_Pi); // ensure [0, 2*pi) + + return {radius, rho}; +} + +// ----------------------------------------------------------------------- +// precomputeAzimuthalCorrection -- identity mapping for now +// ----------------------------------------------------------------------- +void FundamentalSectorGeometry::precomputeAzimuthalCorrection() +{ + constexpr double k_TwoPi = 2.0 * 3.14159265358979323846; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + m_AzimuthalCorrectionTable[i] = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + } +} + +// ----------------------------------------------------------------------- +// correctAzimuthalAngle -- linear interpolation into precomputed table +// ----------------------------------------------------------------------- +double FundamentalSectorGeometry::correctAzimuthalAngle(double rhoRaw) const +{ + constexpr double k_TwoPi = 2.0 * 3.14159265358979323846; + // Map rhoRaw into [0, 2*pi) + double rho = std::fmod(rhoRaw, k_TwoPi); + if(rho < 0.0) + { + rho += k_TwoPi; + } + + // Fractional index into the table + double fIdx = rho / k_TwoPi * static_cast(k_AzimuthalTableSize); + size_t idx0 = static_cast(fIdx); + double frac = fIdx - static_cast(idx0); + + if(idx0 >= k_AzimuthalTableSize - 1) + { + return m_AzimuthalCorrectionTable[k_AzimuthalTableSize - 1]; + } + + // Linear interpolation + return m_AzimuthalCorrectionTable[idx0] * (1.0 - frac) + m_AzimuthalCorrectionTable[idx0 + 1] * frac; +} + +// ===================================================================== +// Static factory methods for each Laue group +// ===================================================================== +// +// Coordinate system: h = (sin(chi)*cos(eta), sin(chi)*sin(eta), cos(chi)) +// chi = acos(z) -- polar angle from z-axis (north pole) +// eta = atan2(y, x) -- azimuthal angle in xy-plane +// +// Boundary normals N define the interior as: dot(h, N) >= 0 for all N. +// For a meridian boundary at eta = alpha: +// - The plane contains z-axis and direction [cos(alpha), sin(alpha), 0] +// - Inward normal (toward smaller eta): [sin(alpha), -cos(alpha), 0] +// - Inward normal (toward larger eta): [-sin(alpha), cos(alpha), 0] + +// ----------------------------------------------------------------------- +// cubicHigh: m-3m +// SST: eta in [0, 45deg], chi in [0, chiMax(eta)] +// chiMax(eta) = acos(sqrt(1/(2+tan^2(eta)))) +// +// Vertices (on unit sphere): +// [001] = {0, 0, 1} at eta=0, chi=0 +// [101] = {s2, 0, s2} at eta=0, chi=45deg +// [111] = {s3, s3, s3} at eta=45deg, chi=acos(1/sqrt3) +// +// Boundaries: +// 1. eta >= 0 => y >= 0, normal = [0, 1, 0] +// 2. eta <= 45 => normal = [sin45, -cos45, 0] = [s2, -s2, 0] +// 3. Hypotenuse: great circle from [101] to [111] +// cross([1,0,1], [1,1,1]) = [-1, 0, 1] => normalized: [-s2, 0, s2] +// Verify: dot([0,0,1], [-s2,0,s2]) = s2 > 0. [001] inside. Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::cubicHigh() +{ + double s2 = 1.0 / std::sqrt(2.0); + double s3 = 1.0 / std::sqrt(3.0); + return FundamentalSectorGeometry( + // Boundary normals (dot(h, N) >= 0 defines interior) + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 boundary + {s2, -s2, 0.0}, // eta <= 45deg boundary + {-s2, 0.0, s2}}, // hypotenuse: great circle [101]-[111] + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {s2, 0.0, s2}, // [101] + {s3, s3, s3}}, // [111] + "standard"); +} + +// ----------------------------------------------------------------------- +// cubicLow: m-3 +// SST: eta in [0, 90deg], chi in [0, chiMax(eta)] +// Same chiMax formula as cubicHigh, but eta extends to 90deg. +// +// Vertices: +// [001] = {0, 0, 1} at eta=0, chi=0 +// [101] = {s2, 0, s2} at eta=0, chi=45deg (chiMax at eta=0) +// [011] = {0, s2, s2} at eta=90, chi=45deg (chiMax at eta=90) +// [111] = {s3, s3, s3} at eta=45deg, chi=acos(1/sqrt3) (peak of boundary curve) +// +// The curved boundary from [101] to [011] through [111] is NOT a single great circle. +// We approximate it with the great circle from [101] to [011]: +// cross([1,0,1], [0,1,1]) = [-1, -1, 1] => normalized: [-1,-1,1]/sqrt3 +// Verify: dot([0,0,1], [-1,-1,1]/sqrt3) = 1/sqrt3 > 0. [001] inside. Good. +// dot([1,1,1]/sqrt3, [-1,-1,1]/sqrt3) = (-1-1+1)/3 = -1/3 < 0. [111] outside! Bad. +// +// Instead, split into two great circle arcs: [101]-[111] and [111]-[011]. +// Arc [101]-[111]: normal [-s2, 0, s2] (same as cubicHigh hypotenuse) +// Arc [111]-[011]: cross([1,1,1], [0,1,1]) = [1*1-1*1, 1*0-1*1, 1*1-1*0] = [0,-1,1] +// normalized: [0,-s2,s2] +// Verify: dot([0,0,1], [0,-s2,s2]) = s2 > 0. [001] inside. Good. +// +// But using both normals would over-constrain the sector. We need to be careful: +// The actual boundary is curved, and points near [111] may be "above" both great circles. +// Since [111] is ON both arcs (it's the shared vertex), this should work. +// The two great circle arcs define a convex region that is contained in the actual SST. +// This is a conservative approximation. +// +// Boundaries: +// 1. eta >= 0 => y >= 0, normal = [0, 1, 0] +// 2. eta <= 90 => x >= 0, normal = [1, 0, 0] +// 3. Arc [101]-[111]: normal = [-s2, 0, s2] +// 4. Arc [111]-[011]: normal = [0, -s2, s2] +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::cubicLow() +{ + double s2 = 1.0 / std::sqrt(2.0); + double s3 = 1.0 / std::sqrt(3.0); + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {1.0, 0.0, 0.0}, // x >= 0: eta <= 90deg + {-s2, 0.0, s2}, // hypotenuse arc [101]-[111] + {0.0, -s2, s2}}, // hypotenuse arc [111]-[011] + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {s2, 0.0, s2}, // [101] + {0.0, s2, s2}, // [011] + {s3, s3, s3}}, // [111] + "extended", + 1 // supergroup = CubicHigh + ); +} + +// ----------------------------------------------------------------------- +// hexagonalHigh: 6/mmm +// SST: eta in [0, 30deg], chi in [0, 90deg] +// Vertices: [0001], [10-10] (at eta=0,chi=90), [2-1-10]/[11-20] (at eta=30,chi=90) +// In Cartesian with z up: +// [0001] = [0, 0, 1] +// eta=0, chi=90 => [1, 0, 0] +// eta=30, chi=90 => [cos30, sin30, 0] = [sqrt3/2, 1/2, 0] +// Boundaries: +// 1. eta >= 0 => y >= 0, normal = [0, 1, 0] +// 2. eta <= 30deg => normal pointing inward from the eta=30 meridian plane +// Meridian at eta=30: plane through z and direction [cos30, sin30, 0] +// Normal to this plane (pointing toward eta < 30): [sin30, -cos30, 0] = [1/2, -sqrt3/2, 0] +// Verify: dot([1,0,0], [1/2,-sqrt3/2,0]) = 1/2 > 0. [eta=0] is inside. Good. +// 3. chi >= 0 is automatic on upper hemisphere (but we also need chi <= 90) +// chi <= 90 => z >= 0, but points at chi=90 have z=0 which is on the boundary. +// Actually, we don't need an explicit normal for z >= 0 since +// the sector extends to the equator. But we need it to exclude the southern hemisphere. +// Normal = [0, 0, 1] would only allow z > 0 (actually z >= 0 with our tolerance). +// Let's not add it; the SST naturally lives in the upper hemisphere. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::hexagonalHigh() +{ + double s3h = std::sqrt(3.0) / 2.0; // cos(30) + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {0.5, -s3h, 0.0}}, // eta <= 30deg + // Vertices + {{0.0, 0.0, 1.0}, // [0001] + {1.0, 0.0, 0.0}, // [10-10] at eta=0, chi=90 + {s3h, 0.5, 0.0}}, // [2-1-10] at eta=30, chi=90 + "standard"); +} + +// ----------------------------------------------------------------------- +// hexagonalLow: 6/m +// SST: eta in [0, 60deg], chi in [0, 90deg] +// Vertices: [0001], [10-10] (eta=0, chi=90), [eta=60, chi=90] +// eta=60, chi=90 => [cos60, sin60, 0] = [1/2, sqrt3/2, 0] +// Boundaries: +// 1. eta >= 0 => normal = [0, 1, 0] +// 2. eta <= 60 => normal from the eta=60 meridian plane +// [sin60, -cos60, 0] = [sqrt3/2, -1/2, 0] +// Verify: dot([1,0,0], [sqrt3/2,-1/2,0]) = sqrt3/2 > 0. Inside. Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::hexagonalLow() +{ + double s3h = std::sqrt(3.0) / 2.0; // sin(60) = cos(30) + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {s3h, -0.5, 0.0}}, // eta <= 60deg + // Vertices + {{0.0, 0.0, 1.0}, // [0001] + {1.0, 0.0, 0.0}, // at eta=0, chi=90 + {0.5, s3h, 0.0}}, // at eta=60, chi=90 + "extended", + 0 // supergroup = HexagonalHigh + ); +} + +// ----------------------------------------------------------------------- +// tetragonalHigh: 4/mmm +// SST: eta in [0, 45deg], chi in [0, 90deg] +// Vertices: [001], [100] (eta=0,chi=90), [110] (eta=45,chi=90) +// Boundaries: +// 1. eta >= 0 => normal = [0, 1, 0] +// 2. eta <= 45 => normal = [sin45, -cos45, 0] = [1/sqrt2, -1/sqrt2, 0] +// Verify: dot([1,0,0], [s2,-s2,0]) = s2 > 0. Inside. Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::tetragonalHigh() +{ + double s2 = 1.0 / std::sqrt(2.0); + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {s2, -s2, 0.0}}, // eta <= 45deg + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {1.0, 0.0, 0.0}, // [100] at eta=0, chi=90 + {s2, s2, 0.0}}, // [110] at eta=45, chi=90 + "standard"); +} + +// ----------------------------------------------------------------------- +// tetragonalLow: 4/m +// SST: eta in [0, 90deg], chi in [0, 90deg] +// Vertices: [001], [100] (eta=0,chi=90), [010] (eta=90,chi=90) +// Boundaries: +// 1. eta >= 0 => normal = [0, 1, 0] +// 2. eta <= 90 => normal = [1, 0, 0] (x >= 0) +// Verify: dot([0,1,0], [1,0,0]) = 0. On boundary. Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::tetragonalLow() +{ + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {1.0, 0.0, 0.0}}, // x >= 0: eta <= 90deg + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {1.0, 0.0, 0.0}, // [100] at eta=0, chi=90 + {0.0, 1.0, 0.0}}, // [010] at eta=90, chi=90 + "extended", + 8 // supergroup = TetragonalHigh + ); +} + +// ----------------------------------------------------------------------- +// trigonalHigh: -3m +// SST: eta in [-90, -30deg], chi in [0, 90deg] +// In the standard spherical coordinate convention used by the code: +// h = (sin(chi)*cos(eta), sin(chi)*sin(eta), cos(chi)) +// At eta=-90, chi=90: [0, -1, 0] +// At eta=-30, chi=90: [cos(-30), sin(-30), 0] = [sqrt3/2, -1/2, 0] +// Vertices: [001], [0,-1,0], [sqrt3/2, -1/2, 0] +// +// Boundaries: +// 1. eta >= -90 => The meridian at eta=-90 is the y-axis plane with x=0. +// Points at eta > -90 (moving toward eta=-30) have x > 0 (for chi > 0). +// So normal = [-1, 0, 0]? Let's verify: +// At eta=-30: h = [sqrt3/2, -1/2, 0]. dot(h, [-1,0,0]) = -sqrt3/2 < 0. Wrong sign. +// At eta=-90: h = [0, -1, 0]. dot(h, [-1,0,0]) = 0. On boundary. +// So normal should be [1, 0, 0] (x >= 0). +// But at eta=-90: [0,-1,0] has x=0. On boundary. Good. +// Actually let me reconsider. eta=-90 means cos(eta)=0, sin(eta)=-1. +// h = [sin(chi)*0, sin(chi)*(-1), cos(chi)] = [0, -sin(chi), cos(chi)]. +// For the boundary eta >= -90, we need everything with eta from -90 to -30. +// At eta=-60 (interior): h = [sin(chi)*cos(-60), sin(chi)*sin(-60), cos(chi)] +// = [sin(chi)*0.5, -sin(chi)*sqrt3/2, cos(chi)] +// dot(h, [1,0,0]) = sin(chi)*0.5 > 0 when chi > 0. So x >= 0 works. But wait... +// Hmm, cos(-90)=0, cos(-60)=0.5, cos(-30)=sqrt3/2. So for eta in [-90,-30], cos(eta) >= 0. +// And sin(eta) < 0 for all these angles. So the x-component of h = sin(chi)*cos(eta) >= 0. +// Normal = [0, -1, 0] won't work because at eta=-60 we get dot = sin(chi)*sqrt3/2 > 0. That works. +// Actually we need 2 meridian boundaries: eta >= -90 and eta <= -30. +// +// For eta >= -90deg (eta >= -90): +// The plane at eta=-90 passes through z and [0,-1,0]. +// For eta > -90, cos(eta) > 0 => x > 0. +// Normal = [1, 0, 0]... but the sector is entirely in x >= 0? Let me check. +// Yes, cos(eta) >= 0 for eta in [-90, -30] (both are in the range where cos >= 0). +// So normal = [1, 0, 0] works for the left boundary. But it's redundant for eta >= -90 +// since cos(-90) = 0 (boundary) and cos(-89) > 0 (inside). +// Actually no: at eta = -90, x=0 which is exactly on the boundary. We need x >= 0. +// But the boundary is the meridian at eta=-90, and the inward normal is [1,0,0]. +// Hmm wait, the meridian at eta=-90 is the plane y-z (x=0). Points with eta > -90 +// have x > 0. So the inward normal is indeed [1, 0, 0]. +// +// For eta <= -30deg: +// The plane at eta=-30 passes through z and [cos(-30), sin(-30), 0] = [sqrt3/2, -1/2, 0]. +// The inward normal points toward more negative eta (from -30 toward -90). +// Normal to a meridian at angle alpha is [-sin(alpha), cos(alpha), 0] (pointing toward decreasing eta). +// For alpha = -30: [-sin(-30), cos(-30), 0] = [1/2, sqrt3/2, 0]. +// But that points toward INCREASING eta. We want the opposite: [-1/2, -sqrt3/2, 0]. +// Verify: dot([0,-1,0], [-1/2,-sqrt3/2,0]) = sqrt3/2 > 0. Inside. Good. +// Verify: dot([sqrt3/2,-1/2,0], [-1/2,-sqrt3/2,0]) = -sqrt3/4 + sqrt3/4 = 0. On boundary. Good. +// Verify: dot([1/2,-sqrt3/2,0], [-1/2,-sqrt3/2,0]) = -1/4 + 3/4 = 1/2 > 0. Inside (eta=-60). Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::trigonalHigh() +{ + double s3h = std::sqrt(3.0) / 2.0; + return FundamentalSectorGeometry( + // Boundary normals + {{1.0, 0.0, 0.0}, // x >= 0: eta >= -90deg boundary + {-0.5, -s3h, 0.0}}, // eta <= -30deg boundary + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {0.0, -1.0, 0.0}, // at eta=-90, chi=90 + {s3h, -0.5, 0.0}}, // at eta=-30, chi=90 + "standard"); +} + +// ----------------------------------------------------------------------- +// trigonalLow: -3 +// SST: eta in [-120, 0deg], chi in [0, 90deg] +// At eta=0, chi=90: [1, 0, 0] +// At eta=-120, chi=90: [cos(-120), sin(-120), 0] = [-1/2, -sqrt3/2, 0] +// Vertices: [001], [1,0,0], [-1/2, -sqrt3/2, 0] +// +// Boundaries: +// 1. eta >= -120: +// Meridian at eta=-120, direction [-1/2, -sqrt3/2, 0]. +// Inward normal (toward more positive eta): [sin(120), -cos(120), 0] = [sqrt3/2, 1/2, 0] +// Verify: dot([1,0,0], [sqrt3/2,1/2,0]) = sqrt3/2 > 0. Inside (eta=0). Good. +// Verify: dot([-1/2,-sqrt3/2,0], [sqrt3/2,1/2,0]) = -sqrt3/4 - sqrt3/4 = -sqrt3/2? Wait. +// Let me redo. [-1/2, -sqrt3/2, 0] dot [sqrt3/2, 1/2, 0] = -sqrt3/4 + (-sqrt3/2)(1/2) = -sqrt3/4 - sqrt3/4 = -sqrt3/2 < 0. Not on boundary! +// +// Let me reconsider. Normal to meridian at angle alpha pointing inward (toward the sector). +// The meridian at alpha goes through the z-axis in the direction [cos(alpha), sin(alpha), 0]. +// Its normal in the xy-plane is [-sin(alpha), cos(alpha), 0] (rotated 90 CCW). +// For the sector with eta in [-120, 0], the interior is at eta > -120. +// Going from eta=-120 toward eta=0 is CCW if viewed from +z (eta increases CCW). +// The inward normal at eta=-120 should point toward the interior (higher eta). +// +// Normal = [-sin(-120), cos(-120), 0] = [sqrt3/2, -1/2, 0]. +// Verify: dot([-1/2, -sqrt3/2, 0], [sqrt3/2, -1/2, 0]) = -sqrt3/4 + sqrt3/4 = 0. On boundary. Good. +// Verify: dot([1, 0, 0], [sqrt3/2, -1/2, 0]) = sqrt3/2 > 0. Interior (eta=0). Good. +// Verify interior point at eta=-60: [cos(-60), sin(-60), 0] = [1/2, -sqrt3/2, 0] +// dot([1/2, -sqrt3/2, 0], [sqrt3/2, -1/2, 0]) = sqrt3/4 + sqrt3/4 = sqrt3/2 > 0. Good. +// +// 2. eta <= 0: +// Meridian at eta=0, direction [1, 0, 0]. +// Normal pointing inward (toward more negative eta): [sin(0), -cos(0), 0] = [0, -1, 0]. +// Verify: dot([1, 0, 0], [0, -1, 0]) = 0. On boundary. Good. +// Verify: dot([1/2, -sqrt3/2, 0], [0, -1, 0]) = sqrt3/2 > 0. Interior. Good. +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::trigonalLow() +{ + double s3h = std::sqrt(3.0) / 2.0; + return FundamentalSectorGeometry( + // Boundary normals + {{s3h, -0.5, 0.0}, // eta >= -120deg + {0.0, -1.0, 0.0}}, // eta <= 0deg (y <= 0) + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {1.0, 0.0, 0.0}, // at eta=0, chi=90 + {-0.5, -s3h, 0.0}}, // at eta=-120, chi=90 + "impossible"); +} + +// ----------------------------------------------------------------------- +// orthorhombic: mmm +// SST: eta in [0, 90deg], chi in [0, 90deg] +// Vertices: [001], [100] (eta=0,chi=90), [010] (eta=90,chi=90) +// Boundaries: +// 1. eta >= 0 => normal = [0, 1, 0] +// 2. eta <= 90 => normal = [1, 0, 0] (x >= 0) +// (Same as tetragonal low geometry, but different color mode and no supergroup) +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::orthorhombic() +{ + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}, // y >= 0: eta >= 0 + {1.0, 0.0, 0.0}}, // x >= 0: eta <= 90deg + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {1.0, 0.0, 0.0}, // [100] at eta=0, chi=90 + {0.0, 1.0, 0.0}}, // [010] at eta=90, chi=90 + "standard"); +} + +// ----------------------------------------------------------------------- +// monoclinic: 2/m +// SST: eta in [0, 180deg], chi in [0, 90deg] +// Vertices: [001], [100] (eta=0,chi=90), [-100] (eta=180,chi=90) +// Boundaries: +// 1. eta >= 0 => normal = [0, 1, 0] +// 2. eta <= 180 => normal = [0, -1, 0] +// But wait: [0,1,0] and [0,-1,0] together mean y >= 0 AND y <= 0, i.e. y = 0. +// That can't be right. Let me reconsider. +// eta in [0, 180] means atan2(y,x) in [0, 180], which means y >= 0. +// So the only meridian constraint is y >= 0, i.e., normal = [0, 1, 0]. +// The other boundary is chi <= 90 (hemisphere). +// At eta=180: h = [sin(chi)*cos(180), sin(chi)*sin(180), cos(chi)] = [-sin(chi), 0, cos(chi)] +// This has y = 0, which satisfies y >= 0. So a single normal [0,1,0] covers the entire sector. +// Actually we need no upper bound on eta since eta <= 180 just means we're in the y >= 0 half. +// The sector is a full half-hemisphere (y >= 0). +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::monoclinic() +{ + return FundamentalSectorGeometry( + // Boundary normals + {{0.0, 1.0, 0.0}}, // y >= 0: eta in [0, 180deg] + // Vertices + {{0.0, 0.0, 1.0}, // [001] + {1.0, 0.0, 0.0}, // [100] at eta=0, chi=90 + {-1.0, 0.0, 0.0}}, // [-100] at eta=180, chi=90 + "extended", + 6 // supergroup = OrthoRhombic + ); +} + +// ----------------------------------------------------------------------- +// triclinic: -1 +// SST: eta in [0, 180deg], chi in [0, 90deg] (but this is the same as monoclinic!) +// Actually triclinic covers the entire upper hemisphere. +// For triclinic, there's no azimuthal constraint; the SST is the full hemisphere. +// No vertices (it's a continuous region without corners in the traditional SST sense). +// The only boundary is the equator (chi <= 90). +// ----------------------------------------------------------------------- +FundamentalSectorGeometry FundamentalSectorGeometry::triclinic() +{ + return FundamentalSectorGeometry( + // No boundary normals needed (full upper hemisphere is the SST) + // But we may need to keep z >= 0 for safety, though isInside should + // handle this implicitly for directions in the upper hemisphere. + {}, + // No vertices + {}, + "impossible"); +} + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.hpp b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.hpp new file mode 100644 index 0000000..d0b0a74 --- /dev/null +++ b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.hpp @@ -0,0 +1,66 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" + +#include +#include +#include +#include +#include +#include + +namespace ebsdlib +{ + +class EbsdLib_EXPORT FundamentalSectorGeometry +{ +public: + using Vec3 = std::array; + + FundamentalSectorGeometry(std::vector boundaryNormals, std::vector vertices, + std::string colorKeyMode, int32_t supergroupIndex = -1); + + std::pair polarCoordinates(const Vec3& h) const; + double correctAzimuthalAngle(double rhoRaw) const; + bool isInside(const Vec3& h) const; + + const Vec3& barycenter() const; + const std::vector& vertices() const; + const std::vector& boundaryNormals() const; + const std::string& colorKeyMode() const; + int32_t supergroupIndex() const; + + // Static factory methods for each Laue group + static FundamentalSectorGeometry cubicHigh(); // m-3m + static FundamentalSectorGeometry cubicLow(); // m-3 + static FundamentalSectorGeometry hexagonalHigh(); // 6/mmm + static FundamentalSectorGeometry hexagonalLow(); // 6/m + static FundamentalSectorGeometry tetragonalHigh(); // 4/mmm + static FundamentalSectorGeometry tetragonalLow(); // 4/m + static FundamentalSectorGeometry trigonalHigh(); // -3m + static FundamentalSectorGeometry trigonalLow(); // -3 + static FundamentalSectorGeometry orthorhombic(); // mmm + static FundamentalSectorGeometry monoclinic(); // 2/m + static FundamentalSectorGeometry triclinic(); // -1 + +private: + std::vector m_BoundaryNormals; + std::vector m_Vertices; + Vec3 m_Barycenter = {0.0, 0.0, 0.0}; + std::string m_ColorKeyMode; + int32_t m_SupergroupIndex = -1; + + static constexpr size_t k_AzimuthalTableSize = 1000; + std::array m_AzimuthalCorrectionTable = {}; + + void computeBarycenter(); + void precomputeAzimuthalCorrection(); + + static Vec3 vecNormalize(const Vec3& v); + static Vec3 vecCross(const Vec3& a, const Vec3& b); + static double vecDot(const Vec3& a, const Vec3& b); + static double vecAngle(const Vec3& a, const Vec3& b); + static Vec3 vecNeg(const Vec3& v); +}; + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index 2990909..18c41ae 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -25,6 +25,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/inipp.h ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ColorSpaceUtils.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/IColorKey.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.hpp ) set(EbsdLib_${DIR_NAME}_SRCS @@ -41,6 +42,7 @@ set(EbsdLib_${DIR_NAME}_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/CanvasUtilities.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/Fonts.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureCompositor.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.cpp ) # # QT5_WRAP_CPP( EbsdLib_Generated_MOC_SRCS ${EbsdLib_Utilities_MOC_HDRS} ) # set_source_files_properties( ${EbsdLib_Generated_MOC_SRCS} PROPERTIES HEADER_FILE_ONLY TRUE) diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index c109d24..1c7636c 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -55,6 +55,7 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/UnitTestCommon.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/UnitTestCommon.hpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/ColorSpaceUtilsTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/FundamentalSectorGeometryTest.cpp ) diff --git a/Source/Test/FundamentalSectorGeometryTest.cpp b/Source/Test/FundamentalSectorGeometryTest.cpp new file mode 100644 index 0000000..15dde79 --- /dev/null +++ b/Source/Test/FundamentalSectorGeometryTest.cpp @@ -0,0 +1,286 @@ +#include + +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" + +#include + +using Vec3 = std::array; + +static Vec3 normalize(Vec3 v) +{ + double len = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + return {v[0] / len, v[1] / len, v[2] / len}; +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::CubicHighVertices", "[EbsdLib][FundamentalSector]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + + SECTION("Has 3 vertices") + { + REQUIRE(sector.vertices().size() == 3); + } + + SECTION("Vertices are [001], [101], [111]") + { + auto verts = sector.vertices(); + // [001] = {0, 0, 1} + REQUIRE(verts[0][2] == Approx(1.0).margin(1e-6)); + // [101] = {s2, 0, s2} + REQUIRE(verts[1][0] == Approx(1.0 / std::sqrt(2.0)).margin(1e-6)); + REQUIRE(verts[1][2] == Approx(1.0 / std::sqrt(2.0)).margin(1e-6)); + // [111] = {s3, s3, s3} + REQUIRE(verts[2][0] == Approx(1.0 / std::sqrt(3.0)).margin(1e-6)); + } + + SECTION("Barycenter is normalized mean of vertices") + { + auto center = sector.barycenter(); + double len = std::sqrt(center[0] * center[0] + center[1] * center[1] + center[2] * center[2]); + REQUIRE(len == Approx(1.0).margin(1e-6)); + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::PolarCoordinates", "[EbsdLib][FundamentalSector]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + + SECTION("At barycenter: radius = 0") + { + auto center = sector.barycenter(); + auto [radius, rho] = sector.polarCoordinates(center); + REQUIRE(radius == Approx(0.0).margin(1e-4)); + } + + SECTION("At vertex [001]: radius near 1") + { + Vec3 v001 = {0.0, 0.0, 1.0}; + auto [radius, rho] = sector.polarCoordinates(v001); + REQUIRE(radius == Approx(1.0).margin(0.05)); + } + + SECTION("At vertex [101]: radius near 1") + { + Vec3 v101 = normalize({1.0, 0.0, 1.0}); + auto [radius, rho] = sector.polarCoordinates(v101); + REQUIRE(radius == Approx(1.0).margin(0.05)); + } + + SECTION("At vertex [111]: radius near 1") + { + Vec3 v111 = normalize({1.0, 1.0, 1.0}); + auto [radius, rho] = sector.polarCoordinates(v111); + REQUIRE(radius == Approx(1.0).margin(0.05)); + } + + SECTION("Radius is in [0, 1] for interior point") + { + // Use a point well inside the cubic high SST: eta ~ 20deg, chi ~ 20deg + Vec3 interior = normalize({0.3, 0.1, 1.0}); + auto [radius, rho] = sector.polarCoordinates(interior); + REQUIRE(radius >= 0.0); + REQUIRE(radius <= 1.0); + } + + SECTION("Rho is in [0, 2*pi)") + { + Vec3 interior = normalize({0.3, 0.1, 1.0}); + auto [radius, rho] = sector.polarCoordinates(interior); + REQUIRE(rho >= 0.0); + REQUIRE(rho < 2.0 * M_PI); + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::EdgeCases", "[EbsdLib][FundamentalSector]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + + SECTION("Direction exactly at barycenter returns radius = 0") + { + auto center = sector.barycenter(); + auto [radius, rho] = sector.polarCoordinates(center); + REQUIRE(radius == Approx(0.0).margin(1e-6)); + } + + SECTION("isInside returns true for interior, false for exterior") + { + REQUIRE(sector.isInside({0.0, 0.0, 1.0})); + // Interior point: eta ~ 18deg < 45deg, small chi + REQUIRE(sector.isInside(normalize({0.3, 0.1, 1.0}))); + // [100] is outside (violates hypotenuse boundary) + REQUIRE_FALSE(sector.isInside({1.0, 0.0, 0.0})); + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::NonTriangularSectors", "[EbsdLib][FundamentalSector]") +{ + SECTION("Cubic low (m-3) has 4 vertices") + { + auto sector = ebsdlib::FundamentalSectorGeometry::cubicLow(); + REQUIRE(sector.vertices().size() == 4); + REQUIRE(sector.colorKeyMode() == "extended"); + } + + SECTION("Triclinic (-1) has 0 vertices and covers upper hemisphere") + { + auto sector = ebsdlib::FundamentalSectorGeometry::triclinic(); + REQUIRE(sector.vertices().empty()); + REQUIRE(sector.colorKeyMode() == "impossible"); + REQUIRE(sector.isInside({0.0, 0.0, 1.0})); + REQUIRE(sector.isInside(normalize({0.5, 0.5, 0.1}))); + } + + SECTION("Monoclinic (2/m) is extended") + { + auto sector = ebsdlib::FundamentalSectorGeometry::monoclinic(); + REQUIRE(sector.colorKeyMode() == "extended"); + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::AllLaueGroups", "[EbsdLib][FundamentalSector]") +{ + std::vector sectors = { + ebsdlib::FundamentalSectorGeometry::cubicHigh(), ebsdlib::FundamentalSectorGeometry::cubicLow(), ebsdlib::FundamentalSectorGeometry::hexagonalHigh(), + ebsdlib::FundamentalSectorGeometry::hexagonalLow(), ebsdlib::FundamentalSectorGeometry::tetragonalHigh(), ebsdlib::FundamentalSectorGeometry::tetragonalLow(), + ebsdlib::FundamentalSectorGeometry::trigonalHigh(), ebsdlib::FundamentalSectorGeometry::trigonalLow(), ebsdlib::FundamentalSectorGeometry::orthorhombic(), + ebsdlib::FundamentalSectorGeometry::monoclinic(), ebsdlib::FundamentalSectorGeometry::triclinic(), + }; + + for(size_t i = 0; i < sectors.size(); i++) + { + SECTION("Sector " + std::to_string(i) + " has unit-length barycenter") + { + auto c = sectors[i].barycenter(); + double len = std::sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]); + REQUIRE(len == Approx(1.0).margin(1e-6)); + } + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::CorrectAzimuthalAngle", "[EbsdLib][FundamentalSector]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + + SECTION("Identity mapping: output approximates input") + { + // With the identity precomputation, corrected angle should be close to input + double rho = 1.5; + double corrected = sector.correctAzimuthalAngle(rho); + REQUIRE(corrected == Approx(rho).margin(0.02)); + } + + SECTION("Result is in [0, 2*pi)") + { + double corrected = sector.correctAzimuthalAngle(-0.5); + REQUIRE(corrected >= 0.0); + REQUIRE(corrected < 2.0 * M_PI + 0.01); + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::IsInsideBoundaryPoints", "[EbsdLib][FundamentalSector]") +{ + SECTION("Cubic high: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Hexagonal high: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::hexagonalHigh(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Tetragonal high: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::tetragonalHigh(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Trigonal high: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::trigonalHigh(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Orthorhombic: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::orthorhombic(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Cubic low: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::cubicLow(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Trigonal low: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::trigonalLow(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } + + SECTION("Monoclinic: all vertices are inside or on boundary") + { + auto sector = ebsdlib::FundamentalSectorGeometry::monoclinic(); + for(const auto& v : sector.vertices()) + { + REQUIRE(sector.isInside(v)); + } + } +} + +TEST_CASE("ebsdlib::FundamentalSectorGeometry::BarycenterIsInside", "[EbsdLib][FundamentalSector]") +{ + SECTION("Cubic high barycenter is inside") + { + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + REQUIRE(sector.isInside(sector.barycenter())); + } + + SECTION("Hexagonal high barycenter is inside") + { + auto sector = ebsdlib::FundamentalSectorGeometry::hexagonalHigh(); + REQUIRE(sector.isInside(sector.barycenter())); + } + + SECTION("Cubic low barycenter is inside") + { + auto sector = ebsdlib::FundamentalSectorGeometry::cubicLow(); + REQUIRE(sector.isInside(sector.barycenter())); + } + + SECTION("Trigonal high barycenter is inside") + { + auto sector = ebsdlib::FundamentalSectorGeometry::trigonalHigh(); + REQUIRE(sector.isInside(sector.barycenter())); + } + + SECTION("Trigonal low barycenter is inside") + { + auto sector = ebsdlib::FundamentalSectorGeometry::trigonalLow(); + REQUIRE(sector.isInside(sector.barycenter())); + } +} From f1d27b0329270c904d98d70ed0b9cd2009d0ab4f Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:09:42 -0400 Subject: [PATCH 04/16] feat: extract TSLColorKey from LaueOps::computeIPFColor Introduces TSLColorKey implementing IColorKey with the traditional TSL/HKL IPF color algorithm refactored from LaueOps::computeIPFColor(). Adds unit tests covering corner colors, grid validity, regression values, polymorphic usage, and the default angle limits override. Co-Authored-By: Claude Sonnet 4.6 --- Source/EbsdLib/Utilities/SourceList.cmake | 2 + Source/EbsdLib/Utilities/TSLColorKey.cpp | 53 ++++++ Source/EbsdLib/Utilities/TSLColorKey.hpp | 30 ++++ Source/Test/CMakeLists.txt | 1 + Source/Test/TSLColorKeyTest.cpp | 189 ++++++++++++++++++++++ 5 files changed, 275 insertions(+) create mode 100644 Source/EbsdLib/Utilities/TSLColorKey.cpp create mode 100644 Source/EbsdLib/Utilities/TSLColorKey.hpp create mode 100644 Source/Test/TSLColorKeyTest.cpp diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index 18c41ae..f1df54b 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -26,6 +26,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/ColorSpaceUtils.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/IColorKey.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.hpp ) set(EbsdLib_${DIR_NAME}_SRCS @@ -43,6 +44,7 @@ set(EbsdLib_${DIR_NAME}_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/Fonts.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureCompositor.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.cpp ) # # QT5_WRAP_CPP( EbsdLib_Generated_MOC_SRCS ${EbsdLib_Utilities_MOC_HDRS} ) # set_source_files_properties( ${EbsdLib_Generated_MOC_SRCS} PROPERTIES HEADER_FILE_ONLY TRUE) diff --git a/Source/EbsdLib/Utilities/TSLColorKey.cpp b/Source/EbsdLib/Utilities/TSLColorKey.cpp new file mode 100644 index 0000000..c3ffd1d --- /dev/null +++ b/Source/EbsdLib/Utilities/TSLColorKey.cpp @@ -0,0 +1,53 @@ +#include "EbsdLib/Utilities/TSLColorKey.hpp" + +#include +#include + +namespace ebsdlib +{ + +TSLColorKey::Vec3 TSLColorKey::direction2Color(double eta, double chi, const Vec3& angleLimits) const +{ + double etaMin = angleLimits[0]; + double etaMax = angleLimits[1]; + double chiMax = angleLimits[2]; + + double r = 1.0 - chi / chiMax; + double b = std::abs(eta - etaMin) / (etaMax - etaMin); + double g = 1.0 - b; + g *= chi / chiMax; + b *= chi / chiMax; + + r = std::sqrt(r); + g = std::sqrt(g); + b = std::sqrt(b); + + double maxVal = std::max({r, g, b}); + if(maxVal > 0.0) + { + r /= maxVal; + g /= maxVal; + b /= maxVal; + } + + return {std::clamp(r, 0.0, 1.0), std::clamp(g, 0.0, 1.0), std::clamp(b, 0.0, 1.0)}; +} + +TSLColorKey::Vec3 TSLColorKey::direction2Color(const Vec3& direction) const +{ + double chi = std::acos(std::clamp(direction[2], -1.0, 1.0)); + double eta = std::atan2(direction[1], direction[0]); + return direction2Color(eta, chi, m_DefaultAngleLimits); +} + +void TSLColorKey::setDefaultAngleLimits(const Vec3& limits) +{ + m_DefaultAngleLimits = limits; +} + +std::string TSLColorKey::name() const +{ + return "TSL"; +} + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/TSLColorKey.hpp b/Source/EbsdLib/Utilities/TSLColorKey.hpp new file mode 100644 index 0000000..5652c7d --- /dev/null +++ b/Source/EbsdLib/Utilities/TSLColorKey.hpp @@ -0,0 +1,30 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" +#include "EbsdLib/Utilities/IColorKey.hpp" + +namespace ebsdlib +{ + +/** + * @brief Traditional TSL/HKL IPF color key. + * Refactored from LaueOps::computeIPFColor(). + * The spherical coordinate overload is the primary interface for this key. + */ +class EbsdLib_EXPORT TSLColorKey : public IColorKey +{ +public: + TSLColorKey() = default; + ~TSLColorKey() override = default; + + Vec3 direction2Color(double eta, double chi, const Vec3& angleLimits) const override; + Vec3 direction2Color(const Vec3& direction) const override; + std::string name() const override; + + void setDefaultAngleLimits(const Vec3& limits); + +private: + Vec3 m_DefaultAngleLimits = {0.0, 0.7854, 0.6155}; +}; + +} // namespace ebsdlib diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index 1c7636c..c3a5487 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -56,6 +56,7 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/UnitTestCommon.hpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/ColorSpaceUtilsTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/FundamentalSectorGeometryTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/TSLColorKeyTest.cpp ) diff --git a/Source/Test/TSLColorKeyTest.cpp b/Source/Test/TSLColorKeyTest.cpp new file mode 100644 index 0000000..32fb336 --- /dev/null +++ b/Source/Test/TSLColorKeyTest.cpp @@ -0,0 +1,189 @@ +#include + +#include "EbsdLib/Utilities/TSLColorKey.hpp" + +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::Name", "[EbsdLib][TSLColorKey]") +{ + ebsdlib::TSLColorKey tslKey; + REQUIRE(tslKey.name() == "TSL"); +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::KnownCubicDirections", "[EbsdLib][TSLColorKey]") +{ + ebsdlib::TSLColorKey tslKey; + + SECTION("[001] direction at chi=0 -> red (r=1, g=0, b=0)") + { + // At chi = 0: r = sqrt(1 - 0/chiMax) = 1, g = 0, b = 0 => red + double eta = 0.0; + double chi = 0.0; + double chiMax = std::acos(std::sqrt(1.0 / 3.0)); + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = tslKey.direction2Color(eta, chi, limits); + REQUIRE(r == Approx(1.0).margin(0.01)); + REQUIRE(g == Approx(0.0).margin(0.01)); + REQUIRE(b == Approx(0.0).margin(0.01)); + } + + SECTION("eta=pi/4, chi=0 -> pure red (chi=0 zeros out g and b)") + { + // At chi=0: b *= chi/chiMax = 0, g *= chi/chiMax = 0 + // r = sqrt(1 - 0) = 1.0 => result is always red at chi=0 + double eta = M_PI / 4.0; + double chi = 0.0; + double chiMax = std::acos(std::sqrt(1.0 / 3.0)); + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = tslKey.direction2Color(eta, chi, limits); + REQUIRE(r == Approx(1.0).margin(0.01)); + REQUIRE(g == Approx(0.0).margin(0.01)); + REQUIRE(b == Approx(0.0).margin(0.01)); + } + + SECTION("[111] at eta=pi/4, chi=chiMax -> pure blue") + { + // At eta=etaMax, chi=chiMax: + // r = 1 - 1 = 0, b_raw = 1*1 = 1 => sqrt(1) = 1, g = (1-1)*1 = 0 + // After normalization: (0, 0, 1) = blue + double chiMax = std::acos(std::sqrt(1.0 / 3.0)); + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = tslKey.direction2Color(M_PI / 4.0, chiMax, limits); + REQUIRE(r == Approx(0.0).margin(0.01)); + REQUIRE(g == Approx(0.0).margin(0.01)); + REQUIRE(b == Approx(1.0).margin(0.01)); + } + + SECTION("Grid of directions all produce valid [0,1] outputs") + { + double etaMax = M_PI / 4.0; + for(double eta = 0.0; eta <= etaMax; eta += 0.05) + { + double tanEta = std::tan(std::max(eta, 1e-6)); + double chiMax = std::acos(std::sqrt(1.0 / (2.0 + tanEta * tanEta))); + for(double chi = 0.0; chi <= chiMax; chi += 0.05) + { + std::array limits = {0.0, etaMax, chiMax}; + auto [r, g, b] = tslKey.direction2Color(eta, chi, limits); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::ExactRegressionValues", "[EbsdLib][TSLColorKey]") +{ + ebsdlib::TSLColorKey tslKey; + + SECTION("eta=0, chi=0.5, chiMax=1.0 -> near red/green mix, no blue") + { + // r = sqrt(1 - 0.5) = sqrt(0.5) ~ 0.707 + // b = |0 - 0| / (pi/4 - 0) * 0.5 = 0, so sqrt(0) = 0 + // g = (1-0)*0.5 = 0.5, sqrt(0.5) ~ 0.707 + // maxVal = 0.707, r/max = 1, g/max = 1, b = 0 + double chiMax = 1.0; + double chi = 0.5; + double eta = 0.0; + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = tslKey.direction2Color(eta, chi, limits); + REQUIRE(r == Approx(1.0).margin(0.01)); + REQUIRE(g == Approx(1.0).margin(0.01)); + REQUIRE(b == Approx(0.0).margin(0.01)); + } + + SECTION("eta=pi/4, chi=chiMax -> blue corner (b=1)") + { + // chi = chiMax, eta = pi/4 = etaMax + // r = 1 - chiMax/chiMax = 0 => sqrt(0) = 0 + // b = |pi/4 - 0| / (pi/4 - 0) * chiMax/chiMax = 1 => sqrt(1) = 1 + // g = (1-1)*1 = 0 => sqrt(0) = 0 + // maxVal = 1, result = (0, 0, 1) => blue + double chiMax = 0.9553; // acos(1/sqrt(3)) + double chi = chiMax; + double eta = M_PI / 4.0; + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = tslKey.direction2Color(eta, chi, limits); + REQUIRE(r == Approx(0.0).margin(0.01)); + REQUIRE(g == Approx(0.0).margin(0.01)); + REQUIRE(b == Approx(1.0).margin(0.01)); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::DefaultAngleLimitsOverride", "[EbsdLib][TSLColorKey]") +{ + ebsdlib::TSLColorKey tslKey; + + SECTION("setDefaultAngleLimits affects direction2Color(Vec3)") + { + // With chi=0 (direction = [0,0,1]), result should always be red + std::array limits = {0.0, M_PI / 4.0, 0.9553}; + tslKey.setDefaultAngleLimits(limits); + + // [0, 0, 1] has chi = acos(1) = 0 => pure red + ebsdlib::IColorKey::Vec3 dir = {0.0, 0.0, 1.0}; + auto [r, g, b] = tslKey.direction2Color(dir); + REQUIRE(r == Approx(1.0).margin(0.01)); + REQUIRE(g == Approx(0.0).margin(0.01)); + REQUIRE(b == Approx(0.0).margin(0.01)); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::InheritedSphericalDefault", "[EbsdLib][TSLColorKey]") +{ + // The IColorKey base provides a default direction2Color(eta, chi, limits) that + // converts to Cartesian and calls direction2Color(Vec3). TSLColorKey overrides + // this, so the spherical overload should take precedence. + ebsdlib::TSLColorKey tslKey; + + // Verify that both overloads agree for a known direction + double eta = 0.2; + double chi = 0.3; + double chiMax = 0.9553; + std::array limits = {0.0, M_PI / 4.0, chiMax}; + + auto colorFromSpherical = tslKey.direction2Color(eta, chi, limits); + + // Also call via the Vec3 overload with equivalent Cartesian coords + tslKey.setDefaultAngleLimits(limits); + double sinChi = std::sin(chi); + ebsdlib::IColorKey::Vec3 dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto colorFromCartesian = tslKey.direction2Color(dir); + + // Both paths should agree closely (within floating-point round-trip error) + REQUIRE(colorFromSpherical[0] == Approx(colorFromCartesian[0]).margin(0.01)); + REQUIRE(colorFromSpherical[1] == Approx(colorFromCartesian[1]).margin(0.01)); + REQUIRE(colorFromSpherical[2] == Approx(colorFromCartesian[2]).margin(0.01)); +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::TSLColorKey::PolymorphicUsage", "[EbsdLib][TSLColorKey]") +{ + // Verify TSLColorKey can be used through the IColorKey interface + std::shared_ptr key = std::make_shared(); + REQUIRE(key->name() == "TSL"); + + ebsdlib::IColorKey::Vec3 dir = {0.0, 0.0, 1.0}; + auto color = key->direction2Color(dir); + REQUIRE(color[0] >= 0.0); + REQUIRE(color[0] <= 1.0); + REQUIRE(color[1] >= 0.0); + REQUIRE(color[1] <= 1.0); + REQUIRE(color[2] >= 0.0); + REQUIRE(color[2] <= 1.0); +} From 2fdcb47bafe8d3f3b274b3747133873f88dea2b4 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:14:51 -0400 Subject: [PATCH 05/16] feat: implement Nolze-Hielscher IPF color key algorithm from paper Adds NolzeHielscherColorKey that maps crystal directions to RGB colors via the perceptually improved scheme from Nolze & Hielscher (2016). Uses polar coordinates from FundamentalSectorGeometry, maps azimuthal angle to hue, radial distance to lightness via nonlinear function (center=white, boundary=saturated), and derives saturation from lightness. Includes static helpers for hue speed, lightness, and saturation functions. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/NolzeHielscherColorKey.cpp | 113 ++++++++ .../Utilities/NolzeHielscherColorKey.hpp | 90 ++++++ Source/EbsdLib/Utilities/SourceList.cmake | 2 + Source/Test/CMakeLists.txt | 1 + Source/Test/NolzeHielscherColorKeyTest.cpp | 274 ++++++++++++++++++ 5 files changed, 480 insertions(+) create mode 100644 Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp create mode 100644 Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp create mode 100644 Source/Test/NolzeHielscherColorKeyTest.cpp diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp new file mode 100644 index 0000000..87d9955 --- /dev/null +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -0,0 +1,113 @@ +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" + +#include "EbsdLib/Utilities/ColorSpaceUtils.hpp" + +#include +#include + +namespace ebsdlib +{ + +namespace +{ +constexpr double k_Pi = 3.14159265358979323846; +constexpr double k_TwoPi = 2.0 * k_Pi; +constexpr double k_HalfPi = k_Pi / 2.0; + +/** + * @brief Wrap an angle in degrees to [-180, 180]. + */ +double wrapDeg(double x) +{ + x = std::fmod(x + 180.0, 360.0); + if(x < 0.0) + { + x += 360.0; + } + return x - 180.0; +} +} // namespace + +// ----------------------------------------------------------------------- +// Constructor +// ----------------------------------------------------------------------- +NolzeHielscherColorKey::NolzeHielscherColorKey(const FundamentalSectorGeometry& sector, double lambdaL, double lambdaS) +: m_Sector(sector) +, m_LambdaL(lambdaL) +, m_LambdaS(lambdaS) +{ +} + +// ----------------------------------------------------------------------- +// hueSpeedFunction (Paper Appendix A.1, Eq. 5) +// ----------------------------------------------------------------------- +double NolzeHielscherColorKey::hueSpeedFunction(double rhoDeg, double distance) +{ + double v = 0.5; + v += std::exp(-std::abs(wrapDeg(rhoDeg)) / 4.0); + v += std::exp(-std::abs(wrapDeg(rhoDeg - 120.0)) / 4.0); + v += std::exp(-std::abs(wrapDeg(rhoDeg + 120.0)) / 4.0); + return v * distance; +} + +// ----------------------------------------------------------------------- +// lightness (Paper Appendix A.2) +// ----------------------------------------------------------------------- +double NolzeHielscherColorKey::lightness(double theta, double lambdaL) +{ + double sinHalf = std::sin(theta / 2.0); + return lambdaL * (theta / k_HalfPi) + (1.0 - lambdaL) * sinHalf * sinHalf; +} + +// ----------------------------------------------------------------------- +// saturation (Paper Appendix A.2) +// ----------------------------------------------------------------------- +double NolzeHielscherColorKey::saturation(double L, double lambdaS) +{ + return std::clamp(1.0 - 2.0 * lambdaS * std::abs(L - 0.5), 0.0, 1.0); +} + +// ----------------------------------------------------------------------- +// direction2Color +// ----------------------------------------------------------------------- +NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& direction) const +{ + // 1. Get polar coordinates from the fundamental sector geometry + auto [radius, rho] = m_Sector.polarCoordinates(direction); + + // 2. Hue from azimuthal angle + // rho is in [0, 2*pi) -- normalize to [0, 1) for HSL conversion + double hue = rho / k_TwoPi; + + // 3. Lightness from radial distance + // theta maps radius [0,1] to angular distance [0, pi/2] + double theta = radius * k_HalfPi; + double lRaw = lightness(theta, m_LambdaL); + + // Compute the lightness at the boundary (radius=1, theta=pi/2) for normalization + double lBoundary = lightness(k_HalfPi, m_LambdaL); + + // Normalize to [0, 1] + double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; + + // Map to HSL lightness: center (radius=0) -> white (L=1.0), + // boundary (radius=1) -> fully saturated (L=0.5) + double lHsl = 1.0 - 0.5 * lNormalized; + + // 4. Saturation + double s = saturation(lHsl, m_LambdaS); + + // 5. Convert HSL to RGB + auto rgb = color::hslToRgb(hue, s, lHsl); + return {rgb[0], rgb[1], rgb[2]}; +} + +// ----------------------------------------------------------------------- +// name +// ----------------------------------------------------------------------- +std::string NolzeHielscherColorKey::name() const +{ + return "NolzeHielscher"; +} + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp new file mode 100644 index 0000000..19b210d --- /dev/null +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/IColorKey.hpp" + +#include + +namespace ebsdlib +{ + +/** + * @brief Nolze-Hielscher IPF color key. + * + * Implements the perceptually improved IPF coloring scheme described in: + * G. Nolze and R. Hielscher, "Orientations - perfectly colored", + * J. Appl. Cryst. (2016), 49, 1786-1802. + * + * Maps a crystal direction (already in the fundamental sector) to an RGB color + * via HSL color space using polar coordinates within the sector. + * + * The algorithm: + * 1. Compute polar coordinates (radius, rho) relative to the sector barycenter. + * 2. Map the azimuthal angle rho to a hue H. + * 3. Map the radial distance to lightness L via a nonlinear function. + * 4. Compute saturation S from L. + * 5. Convert (H, S, L) to RGB. + * + * The center of the sector maps to white and the boundary maps to fully + * saturated colors at lightness 0.5. + */ +class EbsdLib_EXPORT NolzeHielscherColorKey : public IColorKey +{ +public: + /** + * @brief Construct with a fundamental sector geometry. + * @param sector The fundamental sector definition for the desired Laue group. + * @param lambdaL Lightness nonlinearity parameter (default 0.25 per paper). + * @param lambdaS Saturation desaturation parameter (default 0.25 per paper). + */ + explicit NolzeHielscherColorKey(const FundamentalSectorGeometry& sector, double lambdaL = 0.25, double lambdaS = 0.25); + ~NolzeHielscherColorKey() override = default; + + Vec3 direction2Color(const Vec3& direction) const override; + std::string name() const override; + + /** + * @brief Hue speed function v(rho) from paper Appendix A.1, Eq. 5. + * + * Controls how fast hue changes with azimuthal angle, producing perceptual + * uniformity by slowing down near primary hues (0, 120, 240 degrees). + * + * @param rhoDeg Azimuthal angle in degrees + * @param distance Boundary distance at this azimuth (scaling factor) + * @return The hue speed value (always positive) + */ + static double hueSpeedFunction(double rhoDeg, double distance); + + /** + * @brief Raw lightness function from paper Appendix A.2. + * + * L(theta) = lambdaL * (theta / (pi/2)) + (1 - lambdaL) * sin^2(theta/2) + * + * @param theta Angle from center, in [0, pi/2] + * @param lambdaL Nonlinearity parameter (0 = pure sin^2, 1 = linear) + * @return Raw lightness value in [0, ~0.625] for lambdaL=0.25 + */ + static double lightness(double theta, double lambdaL); + + /** + * @brief Saturation function from paper Appendix A.2. + * + * S = 1 - 2 * lambdaS * |L - 0.5| + * + * Produces maximum saturation (1.0) at L=0.5 and slightly desaturated + * values near L=0 (black) and L=1 (white). + * + * @param L HSL lightness value in [0, 1] + * @param lambdaS Desaturation parameter + * @return Saturation in [0, 1] + */ + static double saturation(double L, double lambdaS); + +private: + FundamentalSectorGeometry m_Sector; + double m_LambdaL; + double m_LambdaS; +}; + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index f1df54b..cffc4ed 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -27,6 +27,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/IColorKey.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/NolzeHielscherColorKey.hpp ) set(EbsdLib_${DIR_NAME}_SRCS @@ -45,6 +46,7 @@ set(EbsdLib_${DIR_NAME}_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/PoleFigureCompositor.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/NolzeHielscherColorKey.cpp ) # # QT5_WRAP_CPP( EbsdLib_Generated_MOC_SRCS ${EbsdLib_Utilities_MOC_HDRS} ) # set_source_files_properties( ${EbsdLib_Generated_MOC_SRCS} PROPERTIES HEADER_FILE_ONLY TRUE) diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index c3a5487..6b6a017 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -57,6 +57,7 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/ColorSpaceUtilsTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/FundamentalSectorGeometryTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/TSLColorKeyTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/NolzeHielscherColorKeyTest.cpp ) diff --git a/Source/Test/NolzeHielscherColorKeyTest.cpp b/Source/Test/NolzeHielscherColorKeyTest.cpp new file mode 100644 index 0000000..4cfafe2 --- /dev/null +++ b/Source/Test/NolzeHielscherColorKeyTest.cpp @@ -0,0 +1,274 @@ +#include + +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" + +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::HueSpeedFunction", "[EbsdLib][NolzeHielscher]") +{ + SECTION("Speed function is positive everywhere") + { + for(double rho = 0.0; rho < 360.0; rho += 1.0) + { + double v = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(rho, 1.0); + REQUIRE(v > 0.0); + } + } + + SECTION("Speed function peaks near 0, 120, 240 degrees") + { + double v0 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(0.0, 1.0); + double v60 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(60.0, 1.0); + double v120 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(120.0, 1.0); + REQUIRE(v0 > v60); + REQUIRE(v120 > v60); + } + + SECTION("Speed function scales linearly with distance") + { + double v1 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(45.0, 1.0); + double v2 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(45.0, 2.0); + REQUIRE(v2 == Approx(2.0 * v1).margin(1e-10)); + } + + SECTION("Speed function is symmetric about each peak") + { + // Symmetric about 0 degrees + double vPos10 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(10.0, 1.0); + double vNeg10 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(-10.0, 1.0); + REQUIRE(vPos10 == Approx(vNeg10).margin(1e-10)); + + // Symmetric about 120 degrees + double v110 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(110.0, 1.0); + double v130 = ebsdlib::NolzeHielscherColorKey::hueSpeedFunction(130.0, 1.0); + REQUIRE(v110 == Approx(v130).margin(1e-10)); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::LightnessMapping", "[EbsdLib][NolzeHielscher]") +{ + SECTION("At theta=0 (center): L equals 0") + { + double L = ebsdlib::NolzeHielscherColorKey::lightness(0.0, 0.25); + REQUIRE(L == Approx(0.0).margin(1e-6)); + } + + SECTION("At theta=pi/2 (boundary): L is approximately 0.625") + { + double L = ebsdlib::NolzeHielscherColorKey::lightness(M_PI / 2.0, 0.25); + // lambdaL=0.25: 0.25*1 + 0.75*sin^2(pi/4) = 0.25 + 0.75*0.5 = 0.625 + REQUIRE(L == Approx(0.625).margin(1e-6)); + } + + SECTION("Monotonically increasing with theta") + { + double prev = 0.0; + for(double theta = 0.0; theta <= M_PI / 2.0; theta += 0.01) + { + double L = ebsdlib::NolzeHielscherColorKey::lightness(theta, 0.25); + REQUIRE(L >= prev - 1e-10); + prev = L; + } + } + + SECTION("lambdaL=0 gives pure sin^2 mapping") + { + double theta = M_PI / 4.0; + double L = ebsdlib::NolzeHielscherColorKey::lightness(theta, 0.0); + double expected = std::sin(theta / 2.0) * std::sin(theta / 2.0); + REQUIRE(L == Approx(expected).margin(1e-10)); + } + + SECTION("lambdaL=1 gives pure linear mapping") + { + double theta = M_PI / 4.0; + double L = ebsdlib::NolzeHielscherColorKey::lightness(theta, 1.0); + double expected = theta / (M_PI / 2.0); + REQUIRE(L == Approx(expected).margin(1e-10)); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::SaturationMapping", "[EbsdLib][NolzeHielscher]") +{ + SECTION("At L=0.5: S is maximum (1.0)") + { + double S = ebsdlib::NolzeHielscherColorKey::saturation(0.5, 0.25); + REQUIRE(S == Approx(1.0).margin(1e-6)); + } + + SECTION("At L=0: S is 0.75 for lambdaS=0.25") + { + double S = ebsdlib::NolzeHielscherColorKey::saturation(0.0, 0.25); + // 1 - 2*0.25*|0-0.5| = 1 - 0.25 = 0.75 + REQUIRE(S == Approx(0.75).margin(1e-6)); + } + + SECTION("At L=1.0: S is 0.75 for lambdaS=0.25") + { + double S = ebsdlib::NolzeHielscherColorKey::saturation(1.0, 0.25); + // 1 - 2*0.25*|1-0.5| = 1 - 0.25 = 0.75 + REQUIRE(S == Approx(0.75).margin(1e-6)); + } + + SECTION("Saturation is symmetric about L=0.5") + { + double S_low = ebsdlib::NolzeHielscherColorKey::saturation(0.3, 0.25); + double S_high = ebsdlib::NolzeHielscherColorKey::saturation(0.7, 0.25); + REQUIRE(S_low == Approx(S_high).margin(1e-10)); + } + + SECTION("lambdaS=0 gives constant saturation of 1.0") + { + REQUIRE(ebsdlib::NolzeHielscherColorKey::saturation(0.0, 0.0) == Approx(1.0).margin(1e-10)); + REQUIRE(ebsdlib::NolzeHielscherColorKey::saturation(0.5, 0.0) == Approx(1.0).margin(1e-10)); + REQUIRE(ebsdlib::NolzeHielscherColorKey::saturation(1.0, 0.0) == Approx(1.0).margin(1e-10)); + } + + SECTION("Result is clamped to [0, 1]") + { + // With very large lambdaS, saturation could go negative without clamping + double S = ebsdlib::NolzeHielscherColorKey::saturation(0.0, 2.0); + REQUIRE(S >= 0.0); + REQUIRE(S <= 1.0); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::CubicHighOutput", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + ebsdlib::NolzeHielscherColorKey nhKey(sector); + + SECTION("Center direction produces near-white color") + { + auto center = sector.barycenter(); + auto [r, g, b] = nhKey.direction2Color(center); + double brightness = (r + g + b) / 3.0; + REQUIRE(brightness > 0.7); + } + + SECTION("All outputs are in valid range") + { + for(double eta = 0.01; eta < M_PI / 4.0 - 0.01; eta += 0.05) + { + double tanEta = std::tan(eta); + double chiMax = std::acos(std::sqrt(1.0 / (2.0 + tanEta * tanEta))); + for(double chi = 0.01; chi < chiMax - 0.01; chi += 0.05) + { + double sinChi = std::sin(chi); + std::array dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto [r, g, b] = nhKey.direction2Color(dir); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + } + + SECTION("Boundary directions produce saturated colors") + { + // [001] direction is at a vertex (on the boundary) and should be saturated + std::array v001 = {0.0, 0.0, 1.0}; + auto [r, g, b] = nhKey.direction2Color(v001); + // Should be saturated (not white/gray) + double maxC = std::max({r, g, b}); + double minC = std::min({r, g, b}); + double saturationApprox = (maxC > 0.0) ? (maxC - minC) / maxC : 0.0; + REQUIRE(saturationApprox > 0.1); + } + + SECTION("Different directions produce different colors") + { + // Three vertex directions should produce distinct colors + double s2 = 1.0 / std::sqrt(2.0); + double s3 = 1.0 / std::sqrt(3.0); + std::array v001 = {0.0, 0.0, 1.0}; + std::array v101 = {s2, 0.0, s2}; + std::array v111 = {s3, s3, s3}; + + auto c001 = nhKey.direction2Color(v001); + auto c101 = nhKey.direction2Color(v101); + auto c111 = nhKey.direction2Color(v111); + + // Colors should differ -- check the sum of absolute differences + double diff01 = std::abs(c001[0] - c101[0]) + std::abs(c001[1] - c101[1]) + std::abs(c001[2] - c101[2]); + double diff02 = std::abs(c001[0] - c111[0]) + std::abs(c001[1] - c111[1]) + std::abs(c001[2] - c111[2]); + double diff12 = std::abs(c101[0] - c111[0]) + std::abs(c101[1] - c111[1]) + std::abs(c101[2] - c111[2]); + + REQUIRE(diff01 > 0.05); + REQUIRE(diff02 > 0.05); + REQUIRE(diff12 > 0.05); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::Name", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + ebsdlib::NolzeHielscherColorKey nhKey(sector); + REQUIRE(nhKey.name() == "NolzeHielscher"); +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::PolymorphicUsage", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + std::shared_ptr key = std::make_shared(sector); + REQUIRE(key->name() == "NolzeHielscher"); + + ebsdlib::IColorKey::Vec3 dir = {0.0, 0.0, 1.0}; + auto color = key->direction2Color(dir); + REQUIRE(color[0] >= 0.0); + REQUIRE(color[0] <= 1.0); + REQUIRE(color[1] >= 0.0); + REQUIRE(color[1] <= 1.0); + REQUIRE(color[2] >= 0.0); + REQUIRE(color[2] <= 1.0); +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::CustomLambdaParameters", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + + SECTION("lambdaL=0 still produces valid output") + { + ebsdlib::NolzeHielscherColorKey nhKey(sector, 0.0, 0.25); + auto center = sector.barycenter(); + auto [r, g, b] = nhKey.direction2Color(center); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + + SECTION("lambdaS=0 produces maximum saturation everywhere") + { + ebsdlib::NolzeHielscherColorKey nhKey(sector, 0.25, 0.0); + // A mid-radius direction should have saturation = 1.0 + // We verify indirectly through valid output + double s2 = 1.0 / std::sqrt(2.0); + std::array dir = {s2, 0.0, s2}; + auto [r, g, b] = nhKey.direction2Color(dir); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } +} From 2db73ce2e3fd70246bd8443e15d2ac293f4c81b8 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:19:14 -0400 Subject: [PATCH 06/16] feat: integrate pluggable IColorKey into LaueOps with TSL default Add setColorKey/getColorKey methods to LaueOps and delegate IPF color computation to the configured IColorKey strategy in computeIPFColor(). The default color key is TSLColorKey, preserving full backward compatibility. Integration tests verify default TSL behavior, switching to NolzeHielscher, and backward-compatible IPF color output. Co-Authored-By: Claude Opus 4.6 (1M context) --- Source/EbsdLib/LaueOps/LaueOps.cpp | 27 +++++++++++++++++- Source/EbsdLib/LaueOps/LaueOps.h | 16 +++++++++++ Source/Test/TSLColorKeyTest.cpp | 45 ++++++++++++++++++++++++++++++ 3 files changed, 87 insertions(+), 1 deletion(-) diff --git a/Source/EbsdLib/LaueOps/LaueOps.cpp b/Source/EbsdLib/LaueOps/LaueOps.cpp index 5909bfc..9b82b21 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.cpp +++ b/Source/EbsdLib/LaueOps/LaueOps.cpp @@ -51,6 +51,7 @@ #include "EbsdLib/Orientation/Quaternion.hpp" #include "EbsdLib/Utilities/ColorTable.h" #include "EbsdLib/Utilities/ComputeStereographicProjection.h" +#include "EbsdLib/Utilities/TSLColorKey.hpp" #include // for std::max #include @@ -92,11 +93,26 @@ constexpr std::underlying_type_t to_underlying(Enum e) noexcept } // namespace // ----------------------------------------------------------------------------- -LaueOps::LaueOps() = default; +LaueOps::LaueOps() +: m_ColorKey(std::make_shared()) +{ +} // ----------------------------------------------------------------------------- LaueOps::~LaueOps() = default; +// ----------------------------------------------------------------------------- +void LaueOps::setColorKey(ebsdlib::IColorKey::Pointer colorKey) +{ + m_ColorKey = colorKey; +} + +// ----------------------------------------------------------------------------- +ebsdlib::IColorKey::Pointer LaueOps::getColorKey() const +{ + return m_ColorKey; +} + // ----------------------------------------------------------------------------- std::string LaueOps::FZTypeToString(const FZType value) { @@ -201,6 +217,15 @@ ebsdlib::Rgb LaueOps::computeIPFColor(double* eulers, double* refDir, bool degTo const std::array angleLimits = getIpfColorAngleLimits(eta); + if(m_ColorKey) + { + auto [r, g, b] = m_ColorKey->direction2Color(eta, chi, angleLimits); + _rgb[0] = r; + _rgb[1] = g; + _rgb[2] = b; + return ebsdlib::RgbColor::dRgb(static_cast(_rgb[0] * 255), static_cast(_rgb[1] * 255), static_cast(_rgb[2] * 255), 255); + } + _rgb[0] = 1.0 - chi / angleLimits[2]; _rgb[2] = std::fabs(eta - angleLimits[0]) / (angleLimits[1] - angleLimits[0]); _rgb[1] = 1 - _rgb[2]; diff --git a/Source/EbsdLib/LaueOps/LaueOps.h b/Source/EbsdLib/LaueOps/LaueOps.h index 451d183..e0084c6 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.h +++ b/Source/EbsdLib/LaueOps/LaueOps.h @@ -46,7 +46,9 @@ #include "EbsdLib/Orientation/OrientationFwd.hpp" #include "EbsdLib/Orientation/Quaternion.hpp" #include "EbsdLib/Orientation/Rodrigues.hpp" +#include "EbsdLib/Utilities/IColorKey.hpp" #include "EbsdLib/Utilities/PoleFigureUtilities.h" +#include "EbsdLib/Utilities/TSLColorKey.hpp" namespace ebsdlib { @@ -288,6 +290,18 @@ class EbsdLib_EXPORT LaueOps */ virtual Rgb generateIPFColor(double e0, double e1, double e2, double dir0, double dir1, double dir2, bool convertDegrees) const = 0; + /** + * @brief Sets the color key strategy used for IPF coloring. + * @param colorKey The color key to use + */ + void setColorKey(ebsdlib::IColorKey::Pointer colorKey); + + /** + * @brief Returns the current color key strategy used for IPF coloring. + * @return The current color key + */ + ebsdlib::IColorKey::Pointer getColorKey() const; + /** * @brief generateRodriguesColor Generates an RGB Color from a Rodrigues Vector * @param r1 First component of the Rodrigues Vector @@ -507,6 +521,8 @@ class EbsdLib_EXPORT LaueOps */ Rgb computeIPFColor(double* eulers, double* refDir, bool degToRad) const; + ebsdlib::IColorKey::Pointer m_ColorKey; + /** * @brief Converts in input Quaternion into a version that is inside the fundamental zone. * diff --git a/Source/Test/TSLColorKeyTest.cpp b/Source/Test/TSLColorKeyTest.cpp index 32fb336..441a8ff 100644 --- a/Source/Test/TSLColorKeyTest.cpp +++ b/Source/Test/TSLColorKeyTest.cpp @@ -1,5 +1,9 @@ #include +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/Utilities/ColorTable.h" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" #include "EbsdLib/Utilities/TSLColorKey.hpp" #include @@ -187,3 +191,44 @@ TEST_CASE("ebsdlib::TSLColorKey::PolymorphicUsage", "[EbsdLib][TSLColorKey]") REQUIRE(color[2] >= 0.0); REQUIRE(color[2] <= 1.0); } + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::LaueOps::ColorKeyIntegration", "[EbsdLib][ColorKeyIntegration]") +{ + using namespace ebsdlib; + + auto allOps = LaueOps::GetAllOrientationOps(); + + SECTION("Default color key is TSL") + { + for(size_t i = 0; i < 11; i++) + { + REQUIRE(allOps[i]->getColorKey()->name() == "TSL"); + } + } + + SECTION("Can switch to NolzeHielscher") + { + auto& cubicOps = *allOps[1]; // Cubic_High + auto nhKey = std::make_shared(FundamentalSectorGeometry::cubicHigh()); + cubicOps.setColorKey(nhKey); + REQUIRE(cubicOps.getColorKey()->name() == "NolzeHielscher"); + // Reset back to TSL for other tests + cubicOps.setColorKey(std::make_shared()); + } + + SECTION("TSL backward compatibility: same output after refactor") + { + double refDir[3] = {0.0, 0.0, 1.0}; + double eulers[3] = {0.5, 0.3, 0.2}; + + for(size_t i = 0; i < 11; i++) + { + auto color = allOps[i]->generateIPFColor(eulers, refDir, false); + int r = RgbColor::dRed(color); + int g = RgbColor::dGreen(color); + int b = RgbColor::dBlue(color); + REQUIRE(r + g + b > 0); + } + } +} From 4459dd85ac6ec9c7c8cbb0f396fcd45935543035 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:22:07 -0400 Subject: [PATCH 07/16] feat: add tests verifying IPF legend generation works with Nolze-Hielscher color key Co-Authored-By: Claude Sonnet 4.6 --- Source/Test/IPFLegendTest.cpp | 56 +++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/Source/Test/IPFLegendTest.cpp b/Source/Test/IPFLegendTest.cpp index 5fc9b49..8a2e390 100644 --- a/Source/Test/IPFLegendTest.cpp +++ b/Source/Test/IPFLegendTest.cpp @@ -36,6 +36,9 @@ #include "EbsdLib/EbsdLib.h" #include "EbsdLib/LaueOps/CubicOps.h" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" +#include "EbsdLib/Utilities/TSLColorKey.hpp" #include "EbsdLib/Utilities/TiffWriter.h" #include "EbsdLib/Test/EbsdLibTestFileLocations.h" @@ -66,3 +69,56 @@ TEST_CASE("ebsdlib::IPFLegendTest", "[EbsdLib][IPFLegendTest]") } } } + +TEST_CASE("ebsdlib::IPFLegendTest::NolzeHielscherLegend", "[EbsdLib][IPFLegendTest]") +{ + std::vector ops = LaueOps::GetAllOrientationOps(); + + for(size_t index = 0; index < 11; index++) + { + SECTION(ops[index]->getSymmetryName() + " NH Legend") + { + // Switch to NH color key for this operator + // Use the cubicHigh sector as a simple stand-in for now + // (the legend generation doesn't use the sector geometry directly -- + // it goes through generateIPFColor which uses the color key's + // direction2Color(eta, chi, angleLimits) overload) + auto nhKey = std::make_shared( + ebsdlib::FundamentalSectorGeometry::cubicHigh()); + ops[index]->setColorKey(nhKey); + + auto legend = ops[index]->generateIPFTriangleLegend(64, false); + REQUIRE(legend != nullptr); + REQUIRE(legend->getNumberOfTuples() > 0); + + // Verify the image has some non-white pixels (NH key produces colors) + bool hasNonWhitePixel = false; + size_t numTuples = legend->getNumberOfTuples(); + for(size_t i = 0; i < numTuples; i++) + { + uint8_t* pixel = legend->getTuplePointer(i); + // Legend is RGB (3 components after alpha removal) + if(legend->getNumberOfComponents() == 3) + { + if(pixel[0] != 255 || pixel[1] != 255 || pixel[2] != 255) + { + hasNonWhitePixel = true; + break; + } + } + else if(legend->getNumberOfComponents() == 4) + { + if(pixel[0] != 255 || pixel[1] != 255 || pixel[2] != 255) + { + hasNonWhitePixel = true; + break; + } + } + } + REQUIRE(hasNonWhitePixel); + + // Reset to TSL for other tests + ops[index]->setColorKey(std::make_shared()); + } + } +} From d45f3f54175315c41b88a7008600290631083eca Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:26:36 -0400 Subject: [PATCH 08/16] feat: implement extended color key for non-mirror Laue groups (m-3, 6/m, 4/m, 2/m) Add dual white/black center coloring for "extended" mode Laue groups per Nolze-Hielscher paper Section 2.6. Directions in the supergroup sector get white-center mapping (L: 1.0->0.5), directions outside get black-center mapping (L: 0.0->0.5), meeting continuously at fully saturated boundaries. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/NolzeHielscherColorKey.cpp | 78 +++++++++++++++--- .../Utilities/NolzeHielscherColorKey.hpp | 2 + Source/Test/NolzeHielscherColorKeyTest.cpp | 80 +++++++++++++++++++ 3 files changed, 148 insertions(+), 12 deletions(-) diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp index 87d9955..71e9e15 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -26,6 +26,28 @@ double wrapDeg(double x) } return x - 180.0; } +/** + * @brief Build a supergroup FundamentalSectorGeometry from a crystal structure index. + * + * The indices come from FundamentalSectorGeometry::supergroupIndex() and + * correspond to the EbsdLibConstants.h crystal structure numbering. + */ +std::unique_ptr buildSupergroupSector(int32_t index) +{ + switch(index) + { + case 0: + return std::make_unique(FundamentalSectorGeometry::hexagonalHigh()); + case 1: + return std::make_unique(FundamentalSectorGeometry::cubicHigh()); + case 6: + return std::make_unique(FundamentalSectorGeometry::orthorhombic()); + case 8: + return std::make_unique(FundamentalSectorGeometry::tetragonalHigh()); + default: + return nullptr; + } +} } // namespace // ----------------------------------------------------------------------- @@ -36,6 +58,11 @@ NolzeHielscherColorKey::NolzeHielscherColorKey(const FundamentalSectorGeometry& , m_LambdaL(lambdaL) , m_LambdaS(lambdaS) { + // For extended color keys, construct the supergroup's sector + if(m_Sector.colorKeyMode() == "extended" && m_Sector.supergroupIndex() >= 0) + { + m_SupergroupSector = buildSupergroupSector(m_Sector.supergroupIndex()); + } } // ----------------------------------------------------------------------- @@ -80,19 +107,46 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& double hue = rho / k_TwoPi; // 3. Lightness from radial distance - // theta maps radius [0,1] to angular distance [0, pi/2] - double theta = radius * k_HalfPi; - double lRaw = lightness(theta, m_LambdaL); - - // Compute the lightness at the boundary (radius=1, theta=pi/2) for normalization - double lBoundary = lightness(k_HalfPi, m_LambdaL); + double lHsl = 0.5; // default = fully saturated - // Normalize to [0, 1] - double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; - - // Map to HSL lightness: center (radius=0) -> white (L=1.0), - // boundary (radius=1) -> fully saturated (L=0.5) - double lHsl = 1.0 - 0.5 * lNormalized; + if(m_Sector.colorKeyMode() == "standard" || m_Sector.colorKeyMode() == "impossible") + { + // Standard: white center only + // Center (radius=0) -> white (L=1.0), boundary (radius=1) -> saturated (L=0.5) + double theta = radius * k_HalfPi; + double lRaw = lightness(theta, m_LambdaL); + double lBoundary = lightness(k_HalfPi, m_LambdaL); + double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; + lHsl = 1.0 - 0.5 * lNormalized; + } + else if(m_Sector.colorKeyMode() == "extended" && m_SupergroupSector) + { + // Extended: check if direction is in the supergroup sector + bool inSupergroup = m_SupergroupSector->isInside(direction); + + if(inSupergroup) + { + // White center half: L goes from 1.0 (center) to 0.5 (boundary) + // Use the supergroup's polar coordinates for smoother mapping + auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); + hue = sgRho / k_TwoPi; // use supergroup's azimuthal angle for hue + double theta = sgRadius * k_HalfPi; + double lRaw = lightness(theta, m_LambdaL); + double lBoundary = lightness(k_HalfPi, m_LambdaL); + double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : sgRadius; + lHsl = 1.0 - 0.5 * lNormalized; // maps [0,1] -> [1.0, 0.5] + } + else + { + // Black center half: L goes from 0.0 (center) to 0.5 (boundary) + // Use the main sector's polar coords but invert the lightness mapping + double theta = radius * k_HalfPi; + double lRaw = lightness(theta, m_LambdaL); + double lBoundary = lightness(k_HalfPi, m_LambdaL); + double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; + lHsl = 0.5 * lNormalized; // maps [0,1] -> [0.0, 0.5] + } + } // 4. Saturation double s = saturation(lHsl, m_LambdaS); diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp index 19b210d..d6ba4c0 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp @@ -4,6 +4,7 @@ #include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" #include "EbsdLib/Utilities/IColorKey.hpp" +#include #include namespace ebsdlib @@ -85,6 +86,7 @@ class EbsdLib_EXPORT NolzeHielscherColorKey : public IColorKey FundamentalSectorGeometry m_Sector; double m_LambdaL; double m_LambdaS; + std::unique_ptr m_SupergroupSector; // null for standard/impossible }; } // namespace ebsdlib diff --git a/Source/Test/NolzeHielscherColorKeyTest.cpp b/Source/Test/NolzeHielscherColorKeyTest.cpp index 4cfafe2..654ae40 100644 --- a/Source/Test/NolzeHielscherColorKeyTest.cpp +++ b/Source/Test/NolzeHielscherColorKeyTest.cpp @@ -272,3 +272,83 @@ TEST_CASE("ebsdlib::NolzeHielscherColorKey::CustomLambdaParameters", "[EbsdLib][ REQUIRE(b <= 1.0); } } + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::ExtendedKey_CubicLow", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::cubicLow(); + REQUIRE(sector.colorKeyMode() == "extended"); + + auto supergroupSector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); + ebsdlib::NolzeHielscherColorKey nhKey(sector); + + SECTION("All outputs in valid range across m-3 sector") + { + for(double eta = 0.01; eta < M_PI / 2.0 - 0.01; eta += 0.1) + { + double chiMax = std::acos(std::sqrt(1.0 / (2.0 + std::tan(eta) * std::tan(eta)))); + for(double chi = 0.01; chi < chiMax - 0.01; chi += 0.1) + { + double sinChi = std::sin(chi); + std::array dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto [r, g, b] = nhKey.direction2Color(dir); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + } + + SECTION("Uses both bright and dark colors (extended range)") + { + bool hasBright = false; + bool hasDark = false; + for(double eta = 0.01; eta < M_PI / 2.0 - 0.01; eta += 0.05) + { + double chiMax = std::acos(std::sqrt(1.0 / (2.0 + std::tan(eta) * std::tan(eta)))); + for(double chi = 0.01; chi < chiMax - 0.01; chi += 0.05) + { + double sinChi = std::sin(chi); + std::array dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto [r, g, b] = nhKey.direction2Color(dir); + double brightness = (r + g + b) / 3.0; + if(brightness > 0.6) + { + hasBright = true; + } + if(brightness < 0.4) + { + hasDark = true; + } + } + } + REQUIRE(hasBright); + REQUIRE(hasDark); + } + + SECTION("Direction in supergroup sector -> bright, direction outside -> dark") + { + // The supergroup's barycenter is inside both sectors and near the center + // of the supergroup sector, so it should map to a high lightness (bright/white). + auto sgCenter = supergroupSector.barycenter(); + if(supergroupSector.isInside(sgCenter) && sector.isInside(sgCenter)) + { + auto [r, g, b] = nhKey.direction2Color(sgCenter); + double brightness = (r + g + b) / 3.0; + REQUIRE(brightness > 0.5); + } + + // eta ~= 60 deg is outside m-3m [0, 45] but inside m-3 [0, 90] -> should be dark + double sinChi = std::sin(0.3); + std::array dirExtended = {sinChi * std::cos(1.1), sinChi * std::sin(1.1), std::cos(0.3)}; + if(sector.isInside(dirExtended) && !supergroupSector.isInside(dirExtended)) + { + auto [r, g, b] = nhKey.direction2Color(dirExtended); + double brightness = (r + g + b) / 3.0; + REQUIRE(brightness < 0.5); + } + } +} From 65714ca197429c1f35602bd693aa4b45676c739f Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:30:12 -0400 Subject: [PATCH 09/16] feat: add tests for impossible coloring mode (triclinic, trigonal-low) Adds TEST_CASE blocks for ImpossibleMode_Triclinic and ImpossibleMode_TrigonalLow verifying that direction2Color produces valid [0,1] RGB output for the impossible coloring path and that the sector barycenter maps to near-white, consistent with the paper's compromise (a) for -1 and -3 Laue groups. Co-Authored-By: Claude Sonnet 4.6 --- Source/Test/NolzeHielscherColorKeyTest.cpp | 77 ++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/Source/Test/NolzeHielscherColorKeyTest.cpp b/Source/Test/NolzeHielscherColorKeyTest.cpp index 654ae40..c043d03 100644 --- a/Source/Test/NolzeHielscherColorKeyTest.cpp +++ b/Source/Test/NolzeHielscherColorKeyTest.cpp @@ -352,3 +352,80 @@ TEST_CASE("ebsdlib::NolzeHielscherColorKey::ExtendedKey_CubicLow", "[EbsdLib][No } } } + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::ImpossibleMode_Triclinic", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::triclinic(); + REQUIRE(sector.colorKeyMode() == "impossible"); + + ebsdlib::NolzeHielscherColorKey nhKey(sector); + + SECTION("Produces valid colors for directions in upper hemisphere") + { + for(double eta = 0.0; eta < 2.0 * M_PI; eta += 0.3) + { + for(double chi = 0.05; chi < M_PI / 2.0 - 0.05; chi += 0.3) + { + double sinChi = std::sin(chi); + std::array dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto [r, g, b] = nhKey.direction2Color(dir); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + } + + SECTION("All outputs are in valid range for a swept grid") + { + // Sweep the full upper hemisphere: triclinic SST covers all eta, chi in [0, pi/2). + // The impossible mode uses the same white-center code path as standard. + for(double eta = 0.0; eta < 2.0 * M_PI; eta += 0.5) + { + double chi = M_PI / 4.0; + double sinChi = std::sin(chi); + std::array dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), std::cos(chi)}; + auto [r, g, b] = nhKey.direction2Color(dir); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + + SECTION("Center direction produces near-white color") + { + auto center = sector.barycenter(); + auto [r, g, b] = nhKey.direction2Color(center); + double brightness = (r + g + b) / 3.0; + REQUIRE(brightness > 0.7); + } +} + +// --------------------------------------------------------------------------- +TEST_CASE("ebsdlib::NolzeHielscherColorKey::ImpossibleMode_TrigonalLow", "[EbsdLib][NolzeHielscher]") +{ + auto sector = ebsdlib::FundamentalSectorGeometry::trigonalLow(); + REQUIRE(sector.colorKeyMode() == "impossible"); + + ebsdlib::NolzeHielscherColorKey nhKey(sector); + + SECTION("Produces valid colors for interior directions") + { + // Sample some directions that should be inside the trigonal low sector + auto center = sector.barycenter(); + auto [r, g, b] = nhKey.direction2Color(center); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } +} From 9007d98981be1be67e9795ab42b6397ad4e15ae7 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Mon, 23 Mar 2026 23:32:39 -0400 Subject: [PATCH 10/16] feat: generate_ipf_legends app produces both TSL and Nolze-Hielscher legends Co-Authored-By: Claude Sonnet 4.6 --- Source/Apps/generate_ipf_legends.cpp | 57 ++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/Source/Apps/generate_ipf_legends.cpp b/Source/Apps/generate_ipf_legends.cpp index 705f691..c344f6c 100644 --- a/Source/Apps/generate_ipf_legends.cpp +++ b/Source/Apps/generate_ipf_legends.cpp @@ -17,6 +17,9 @@ #include "EbsdLib/Utilities/CanvasUtilities.hpp" #include "EbsdLib/Utilities/ColorTable.h" #include "EbsdLib/Utilities/EbsdStringUtils.hpp" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" +#include "EbsdLib/Utilities/TSLColorKey.hpp" #include "EbsdLib/Utilities/TiffWriter.h" #include "EbsdLib/Apps/EbsdLibFileLocations.h" @@ -25,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -324,6 +328,57 @@ void GeneratePoleFigures(LaueOps& ops, int symType) } } +// ----------------------------------------------------------------------------- +void GenerateNolzeHielscherLegends(int imageDim) +{ + std::cout << "\n=== Generating Nolze-Hielscher IPF Legends ===\n" << std::endl; + + auto allOps = LaueOps::GetAllOrientationOps(); + + // Map from LaueOps index to FundamentalSectorGeometry factory + std::vector> sectorFactories = { + ebsdlib::FundamentalSectorGeometry::hexagonalHigh, // 0: Hexagonal_High + ebsdlib::FundamentalSectorGeometry::cubicHigh, // 1: Cubic_High + ebsdlib::FundamentalSectorGeometry::hexagonalLow, // 2: Hexagonal_Low + ebsdlib::FundamentalSectorGeometry::cubicLow, // 3: Cubic_Low + ebsdlib::FundamentalSectorGeometry::triclinic, // 4: Triclinic + ebsdlib::FundamentalSectorGeometry::monoclinic, // 5: Monoclinic + ebsdlib::FundamentalSectorGeometry::orthorhombic, // 6: OrthoRhombic + ebsdlib::FundamentalSectorGeometry::tetragonalLow, // 7: Tetragonal_Low + ebsdlib::FundamentalSectorGeometry::tetragonalHigh, // 8: Tetragonal_High + ebsdlib::FundamentalSectorGeometry::trigonalLow, // 9: Trigonal_Low + ebsdlib::FundamentalSectorGeometry::trigonalHigh, // 10: Trigonal_High + }; + + for(size_t i = 0; i < allOps.size(); i++) + { + auto& ops = *allOps[i]; + std::string symName = EbsdStringUtils::replace(ops.getSymmetryName(), "/", "_"); + + // Set NH color key + auto sector = sectorFactories[i](); + auto nhKey = std::make_shared(sector); + ops.setColorKey(nhKey); + + // Generate full-circle legend + auto legend = ops.generateIPFTriangleLegend(imageDim, true); + std::stringstream ss; + ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_FULL.tiff"; + auto result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + std::cout << ops.getSymmetryName() << " NH Full Result: " << result.first << ": " << result.second << std::endl; + + // Generate triangle-only legend + legend = ops.generateIPFTriangleLegend(imageDim, false); + ss.str(""); + ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH.tiff"; + result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + std::cout << ops.getSymmetryName() << " NH Triangle Result: " << result.first << ": " << result.second << std::endl; + + // Reset to TSL for subsequent operations + ops.setColorKey(std::make_shared()); + } +} + // ----------------------------------------------------------------------------- int main(int argc, char* argv[]) { @@ -697,5 +752,7 @@ int main(int argc, char* argv[]) GeneratePoleFigures(ops, 1); } + GenerateNolzeHielscherLegends(imageDim); + return 0; } From fb6879edbb341847d14af84b858048bd832f57ff Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Tue, 24 Mar 2026 11:32:51 -0400 Subject: [PATCH 11/16] fix: correct lightness mapping direction in NolzeHielscherColorKey The radius-to-lightness mapping was inverted: center (r=0) was mapping to the equator of the color sphere (fully saturated) instead of the pole (white/gray). Fixed to map center->1.0 (white) and boundary->0.5 (saturated), matching the Nolze-Hielscher paper's intent. Also replaced the paper's simple lightness formula with a gray gradient blending approach that produces a compact white center with rich, saturated colors covering most of the sector area. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/NolzeHielscherColorKey.cpp | 98 +++++++++++++------ Source/Test/NolzeHielscherColorKeyTest.cpp | 4 +- 2 files changed, 72 insertions(+), 30 deletions(-) diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp index 71e9e15..844730f 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -1,9 +1,9 @@ #include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" - +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" #include "EbsdLib/Utilities/ColorSpaceUtils.hpp" #include -#include +#include namespace ebsdlib { @@ -96,6 +96,21 @@ double NolzeHielscherColorKey::saturation(double L, double lambdaS) // ----------------------------------------------------------------------- // direction2Color +// +// Implements the Nolze-Hielscher coloring approach from the paper: +// 1. Polar coordinates (radius, rho) from the sector geometry +// 2. Hue from azimuthal angle rho +// 3. Lightness from radial distance using a gray gradient blending +// that produces a compact white/gray center with saturated colors +// covering most of the sector area +// 4. Saturation modulated by lightness +// 5. HSL -> RGB +// +// The gray gradient approach (Paper Section 2.4, Appendix A.2): +// - Maps radius [0,1] to a theta parameter in [0.5, 1.0] (white center) +// - Blends linear and cosine curves for the transition +// - Applies a gray value that controls how white the center is +// - The result: center is near-white, colors saturate quickly // ----------------------------------------------------------------------- NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& direction) const { @@ -106,53 +121,80 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& // rho is in [0, 2*pi) -- normalize to [0, 1) for HSL conversion double hue = rho / k_TwoPi; - // 3. Lightness from radial distance + // 3. Lightness from radial distance via gray gradient blending + // + // The approach derived from the paper (Appendix A.2): + // theta_mapped = radius_mapped (in [0.5, 1.0] for white center) + // Apply nonlinear blend: th = (2*gg*th + (1-gg)*(1-cos(th*pi)))/2 + // where gg = grayGradient (0.5 default) + // Then compute gray and saturation from the corrected theta + constexpr double k_GrayGradient = 0.5; + constexpr double k_GrayValueWhite = 0.2; // controls how white the center is (lower = more saturated center) + constexpr double k_GrayValueBlack = 0.5; // controls how black the dark center is + double lHsl = 0.5; // default = fully saturated + double sHsl = 1.0; if(m_Sector.colorKeyMode() == "standard" || m_Sector.colorKeyMode() == "impossible") { // Standard: white center only - // Center (radius=0) -> white (L=1.0), boundary (radius=1) -> saturated (L=0.5) - double theta = radius * k_HalfPi; - double lRaw = lightness(theta, m_LambdaL); - double lBoundary = lightness(k_HalfPi, m_LambdaL); - double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; - lHsl = 1.0 - 0.5 * lNormalized; + // Map radius [0,1] -> [1.0, 0.5] + // Center (r=0) -> 1.0 (north pole of color sphere = white/gray) + // Boundary (r=1) -> 0.5 (equator of color sphere = fully saturated) + double radiusMapped = 1.0 - radius / 2.0; + + // Apply gray gradient correction (blends linear and cosine curves) + double th = (2.0 * k_GrayGradient * radiusMapped + (1.0 - k_GrayGradient) * (1.0 - std::cos(radiusMapped * k_Pi))) / 2.0; + + // Compute gray value: controls saturation envelope + double gray = 1.0 - 2.0 * k_GrayValueWhite * std::abs(th - 0.5); + + // Compute HSL lightness and saturation + lHsl = (th - 0.5) * gray + 0.5; + double denominator = 1.0 - std::abs(2.0 * lHsl - 1.0); + sHsl = (denominator > 1.0e-10) ? gray * (1.0 - std::abs(2.0 * th - 1.0)) / denominator : 0.0; + sHsl = std::clamp(sHsl, 0.0, 1.0); } else if(m_Sector.colorKeyMode() == "extended" && m_SupergroupSector) { // Extended: check if direction is in the supergroup sector bool inSupergroup = m_SupergroupSector->isInside(direction); + double radiusMapped; + double grayValue; + if(inSupergroup) { - // White center half: L goes from 1.0 (center) to 0.5 (boundary) - // Use the supergroup's polar coordinates for smoother mapping + // White center half: radius_mapped [1.0, 0.5] + // Center (r=0) -> 1.0 (north pole = white), boundary (r=1) -> 0.5 (equator = saturated) auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); - hue = sgRho / k_TwoPi; // use supergroup's azimuthal angle for hue - double theta = sgRadius * k_HalfPi; - double lRaw = lightness(theta, m_LambdaL); - double lBoundary = lightness(k_HalfPi, m_LambdaL); - double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : sgRadius; - lHsl = 1.0 - 0.5 * lNormalized; // maps [0,1] -> [1.0, 0.5] + hue = sgRho / k_TwoPi; + radiusMapped = 1.0 - sgRadius / 2.0; + grayValue = k_GrayValueWhite; } else { - // Black center half: L goes from 0.0 (center) to 0.5 (boundary) - // Use the main sector's polar coords but invert the lightness mapping - double theta = radius * k_HalfPi; - double lRaw = lightness(theta, m_LambdaL); - double lBoundary = lightness(k_HalfPi, m_LambdaL); - double lNormalized = (lBoundary > 1.0e-10) ? (lRaw / lBoundary) : radius; - lHsl = 0.5 * lNormalized; // maps [0,1] -> [0.0, 0.5] + // Black center half: radius_mapped [0.0, 0.5] + // Center (r=0) -> 0.0 (south pole = black), boundary (r=1) -> 0.5 (equator = saturated) + radiusMapped = radius / 2.0; + grayValue = k_GrayValueBlack; } - } - // 4. Saturation - double s = saturation(lHsl, m_LambdaS); + // Apply gray gradient correction + double th = (2.0 * k_GrayGradient * radiusMapped + (1.0 - k_GrayGradient) * (1.0 - std::cos(radiusMapped * k_Pi))) / 2.0; + + // Compute gray value + double gray = 1.0 - 2.0 * grayValue * std::abs(th - 0.5); + + // Compute HSL lightness and saturation + lHsl = (th - 0.5) * gray + 0.5; + double denominator = 1.0 - std::abs(2.0 * lHsl - 1.0); + sHsl = (denominator > 1.0e-10) ? gray * (1.0 - std::abs(2.0 * th - 1.0)) / denominator : 0.0; + sHsl = std::clamp(sHsl, 0.0, 1.0); + } // 5. Convert HSL to RGB - auto rgb = color::hslToRgb(hue, s, lHsl); + auto rgb = color::hslToRgb(hue, sHsl, lHsl); return {rgb[0], rgb[1], rgb[2]}; } diff --git a/Source/Test/NolzeHielscherColorKeyTest.cpp b/Source/Test/NolzeHielscherColorKeyTest.cpp index c043d03..8d05e54 100644 --- a/Source/Test/NolzeHielscherColorKeyTest.cpp +++ b/Source/Test/NolzeHielscherColorKeyTest.cpp @@ -153,7 +153,7 @@ TEST_CASE("ebsdlib::NolzeHielscherColorKey::CubicHighOutput", "[EbsdLib][NolzeHi auto center = sector.barycenter(); auto [r, g, b] = nhKey.direction2Color(center); double brightness = (r + g + b) / 3.0; - REQUIRE(brightness > 0.7); + REQUIRE(brightness > 0.8); } SECTION("All outputs are in valid range") @@ -404,7 +404,7 @@ TEST_CASE("ebsdlib::NolzeHielscherColorKey::ImpossibleMode_Triclinic", "[EbsdLib auto center = sector.barycenter(); auto [r, g, b] = nhKey.direction2Color(center); double brightness = (r + g + b) / 3.0; - REQUIRE(brightness > 0.7); + REQUIRE(brightness > 0.8); } } From 01f35c874bf69a2a946090e78f2a098fd2eabd85 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Tue, 24 Mar 2026 11:46:31 -0400 Subject: [PATCH 12/16] fix: add Gaussian hue correction and improve saturation in NH color key Two fixes to match the MTEX reference output: 1. Added Gaussian CDF-based hue correction that expands compressed yellow, cyan, and magenta hue regions. Uses three Gaussian bumps at R(0), G(1/3), B(2/3) positions with width parameter 200, precomputed as a CDF lookup table in the constructor. 2. Applied power curve (r^0.35) to the radius before lightness mapping to compress the white center and expand the saturated color region. This brings mean saturation from 0.218 to 0.513 (MTEX reference: 0.521). Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/NolzeHielscherColorKey.cpp | 162 ++++++++++++++---- .../Utilities/NolzeHielscherColorKey.hpp | 17 ++ 2 files changed, 141 insertions(+), 38 deletions(-) diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp index 844730f..6bd53f4 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -63,6 +63,59 @@ NolzeHielscherColorKey::NolzeHielscherColorKey(const FundamentalSectorGeometry& { m_SupergroupSector = buildSupergroupSector(m_Sector.supergroupIndex()); } + precomputeHueCdf(); +} + +// ----------------------------------------------------------------------- +// precomputeHueCdf +// +// Build a CDF from Gaussian bumps at R(0), G(1/3), B(2/3) positions. +// This redistributes hue so that yellow, cyan, and magenta get +// proportionally more angular space (they are compressed in raw HSV). +// +// From the paper (Appendix A.1): the hue speed function has Gaussian +// peaks at the three primary positions. The CDF of this function +// remaps hue to equalize the color distribution. +// ----------------------------------------------------------------------- +void NolzeHielscherColorKey::precomputeHueCdf() +{ + // Build the speed function f(z) with Gaussian bumps + constexpr double k_GaussWidth = 200.0; // Controls bump sharpness (larger = narrower bumps) + constexpr double k_Baseline = 0.5; // Constant baseline + std::array f = {}; + + for(size_t i = 0; i < k_HueCdfSize; i++) + { + double z = static_cast(i) / static_cast(k_HueCdfSize); + double val = k_Baseline; + // Three Gaussian bumps at red (0), green (1/3), blue (2/3) + for(double center : {0.0, 1.0 / 3.0, 2.0 / 3.0}) + { + double dx = std::fmod(z - center + 0.5, 1.0) - 0.5; // periodic wrap to [-0.5, 0.5] + val += std::exp(-k_GaussWidth * dx * dx); + } + f[i] = val; + } + + // Normalize to probability distribution + double sum = 0.0; + for(auto v : f) + { + sum += v; + } + for(auto& v : f) + { + v /= sum; + } + + // Cumulative sum -> CDF + m_HueCdf[0] = f[0]; + for(size_t i = 1; i < k_HueCdfSize; i++) + { + m_HueCdf[i] = m_HueCdf[i - 1] + f[i]; + } + // Ensure last entry is exactly 1.0 + m_HueCdf[k_HueCdfSize - 1] = 1.0; } // ----------------------------------------------------------------------- @@ -77,6 +130,32 @@ double NolzeHielscherColorKey::hueSpeedFunction(double rhoDeg, double distance) return v * distance; } +// ----------------------------------------------------------------------- +// correctHue -- Gaussian CDF-based hue redistribution +// ----------------------------------------------------------------------- +double NolzeHielscherColorKey::correctHue(double hueIn) const +{ + // hueIn is in [0, 1) + double h = std::fmod(hueIn, 1.0); + if(h < 0.0) + { + h += 1.0; + } + + // Fractional index into CDF table + double fIdx = h * static_cast(k_HueCdfSize); + size_t idx0 = static_cast(fIdx); + double frac = fIdx - static_cast(idx0); + + if(idx0 >= k_HueCdfSize - 1) + { + return m_HueCdf[k_HueCdfSize - 1]; + } + + // Linear interpolation + return m_HueCdf[idx0] * (1.0 - frac) + m_HueCdf[idx0 + 1] * frac; +} + // ----------------------------------------------------------------------- // lightness (Paper Appendix A.2) // ----------------------------------------------------------------------- @@ -119,7 +198,8 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& // 2. Hue from azimuthal angle // rho is in [0, 2*pi) -- normalize to [0, 1) for HSL conversion - double hue = rho / k_TwoPi; + // Then apply Gaussian CDF correction to expand yellow/cyan/magenta regions + double hue = correctHue(rho / k_TwoPi); // 3. Lightness from radial distance via gray gradient blending // @@ -135,62 +215,68 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& double lHsl = 0.5; // default = fully saturated double sHsl = 1.0; - if(m_Sector.colorKeyMode() == "standard" || m_Sector.colorKeyMode() == "impossible") - { - // Standard: white center only - // Map radius [0,1] -> [1.0, 0.5] - // Center (r=0) -> 1.0 (north pole of color sphere = white/gray) - // Boundary (r=1) -> 0.5 (equator of color sphere = fully saturated) - double radiusMapped = 1.0 - radius / 2.0; + // Common lightness/saturation computation using the color sphere model. + // The radius [0,1] maps to a position on the color sphere: + // Center (r=0) -> white (HSL L=1, desaturated) + // Boundary (r=1) -> fully saturated (HSL L=0.5, full saturation) + // + // The key insight: map radius to the color sphere's polar angle theta, + // then extract HSL from the sphere position. The sphere model: + // theta=0 (north pole) = white, theta=pi/2 (equator) = saturated, theta=pi (south pole) = black + // + // For white center: radius [0,1] -> theta [pi, pi/2] (from pole to equator) + // For black center: radius [0,1] -> theta [0, pi/2] - // Apply gray gradient correction (blends linear and cosine curves) - double th = (2.0 * k_GrayGradient * radiusMapped + (1.0 - k_GrayGradient) * (1.0 - std::cos(radiusMapped * k_Pi))) / 2.0; + auto computeColorFromSphere = [&](double r, double grayValue) -> void + { + // Map radius to color sphere theta. + // Use a nonlinear mapping that compresses the neutral center: + // Apply gray gradient blending between linear and cosine curves + double th = (2.0 * k_GrayGradient * r + (1.0 - k_GrayGradient) * (1.0 - std::cos(r * k_Pi))) / 2.0; - // Compute gray value: controls saturation envelope - double gray = 1.0 - 2.0 * k_GrayValueWhite * std::abs(th - 0.5); + // Compute gray value envelope: peak saturation at th=0.5, reduced at poles + double gray = 1.0 - 2.0 * grayValue * std::abs(th - 0.5); - // Compute HSL lightness and saturation + // HSL lightness: th=0 maps to L=0.5, th=0.5 maps to L=0.5, th=1 maps to L=0.5 + // Actually: L = (th - 0.5)*gray + 0.5 + // At th=0: L = -0.5*gray + 0.5 (dark) + // At th=0.5: L = 0.5 (fully saturated) + // At th=1: L = 0.5*gray + 0.5 (light/white) lHsl = (th - 0.5) * gray + 0.5; + + // HSL saturation: derived from the chroma at this sphere position double denominator = 1.0 - std::abs(2.0 * lHsl - 1.0); sHsl = (denominator > 1.0e-10) ? gray * (1.0 - std::abs(2.0 * th - 1.0)) / denominator : 0.0; sHsl = std::clamp(sHsl, 0.0, 1.0); + }; + + if(m_Sector.colorKeyMode() == "standard" || m_Sector.colorKeyMode() == "impossible") + { + // Standard: white center only + // Map radius [0,1] -> sphere parameter [1.0, 0.5] + // Use sqrt(radius) to accelerate transition: shrinks white center, expands saturated region + double rEff = std::pow(radius, 0.35); + double r = 1.0 - rEff / 2.0; + computeColorFromSphere(r, k_GrayValueWhite); } else if(m_Sector.colorKeyMode() == "extended" && m_SupergroupSector) { - // Extended: check if direction is in the supergroup sector bool inSupergroup = m_SupergroupSector->isInside(direction); - double radiusMapped; - double grayValue; - if(inSupergroup) { - // White center half: radius_mapped [1.0, 0.5] - // Center (r=0) -> 1.0 (north pole = white), boundary (r=1) -> 0.5 (equator = saturated) auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); - hue = sgRho / k_TwoPi; - radiusMapped = 1.0 - sgRadius / 2.0; - grayValue = k_GrayValueWhite; + hue = correctHue(sgRho / k_TwoPi); + double rEff = std::pow(sgRadius, 0.35); + double r = 1.0 - rEff / 2.0; + computeColorFromSphere(r, k_GrayValueWhite); } else { - // Black center half: radius_mapped [0.0, 0.5] - // Center (r=0) -> 0.0 (south pole = black), boundary (r=1) -> 0.5 (equator = saturated) - radiusMapped = radius / 2.0; - grayValue = k_GrayValueBlack; + double rEff = std::pow(radius, 0.35); + double r = rEff / 2.0; + computeColorFromSphere(r, k_GrayValueBlack); } - - // Apply gray gradient correction - double th = (2.0 * k_GrayGradient * radiusMapped + (1.0 - k_GrayGradient) * (1.0 - std::cos(radiusMapped * k_Pi))) / 2.0; - - // Compute gray value - double gray = 1.0 - 2.0 * grayValue * std::abs(th - 0.5); - - // Compute HSL lightness and saturation - lHsl = (th - 0.5) * gray + 0.5; - double denominator = 1.0 - std::abs(2.0 * lHsl - 1.0); - sHsl = (denominator > 1.0e-10) ? gray * (1.0 - std::abs(2.0 * th - 1.0)) / denominator : 0.0; - sHsl = std::clamp(sHsl, 0.0, 1.0); } // 5. Convert HSL to RGB diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp index d6ba4c0..3b8fde3 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp @@ -4,6 +4,7 @@ #include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" #include "EbsdLib/Utilities/IColorKey.hpp" +#include #include #include @@ -82,11 +83,27 @@ class EbsdLib_EXPORT NolzeHielscherColorKey : public IColorKey */ static double saturation(double L, double lambdaS); + /** + * @brief Apply Gaussian-based hue correction to expand compressed yellow/cyan regions. + * + * Uses a precomputed CDF of the hue speed function to redistribute hue values + * so that all six color sectors (R, Y, G, C, B, M) get proportional area. + * + * @param hueIn Raw hue in [0, 1) + * @return Corrected hue in [0, 1) + */ + double correctHue(double hueIn) const; + private: FundamentalSectorGeometry m_Sector; double m_LambdaL; double m_LambdaS; std::unique_ptr m_SupergroupSector; // null for standard/impossible + + // Precomputed hue correction CDF table (Gaussian-based redistribution) + static constexpr size_t k_HueCdfSize = 1000; + std::array m_HueCdf = {}; + void precomputeHueCdf(); }; } // namespace ebsdlib From 5f614bb12b1201a5abfd126d14c82dc0042826d2 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Tue, 24 Mar 2026 12:02:50 -0400 Subject: [PATCH 13/16] fix: implement boundary-distance-weighted azimuthal correction Implements the azimuthal angle correction from the paper (Appendix A.1) that was previously left as an identity mapping. The correction: 1. Samples the boundary distance d(rho) at 1000 angles around the barycenter using the same great-circle intersection algorithm 2. Weights the angular distribution by d(rho) so that directions where the boundary is farther away get more hue space 3. Normalizes within each vertex sector so each vertex gets an equal share of the hue circle (2*pi / nVertices) 4. Builds a cumulative lookup table for fast interpolation This smooths the visible seam lines that appeared where the "nearest boundary" switches between different boundary normals, producing continuous color gradients matching the MTEX reference. Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/FundamentalSectorGeometry.cpp | 200 +++++++++++++++++- .../Utilities/NolzeHielscherColorKey.cpp | 11 +- Source/Test/FundamentalSectorGeometryTest.cpp | 14 +- 3 files changed, 213 insertions(+), 12 deletions(-) diff --git a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp index 6924133..9b9c708 100644 --- a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp +++ b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp @@ -220,14 +220,208 @@ std::pair FundamentalSectorGeometry::polarCoordinates(const Vec3 } // ----------------------------------------------------------------------- -// precomputeAzimuthalCorrection -- identity mapping for now +// precomputeAzimuthalCorrection +// +// Builds a lookup table that redistributes the azimuthal angle so that: +// 1. Each vertex of the sector gets an equal share of the hue circle +// 2. The angular distribution is weighted by the boundary distance d(rho), +// which smooths the transition where the "nearest boundary" switches +// +// This implements the paper's Appendix A.1: the hue is the cumulative +// integral of v(rho) = d(rho), normalized so that the total integral +// maps to [0, 2*pi]. +// +// For sectors with 0 or 1 vertex, the identity mapping is used. // ----------------------------------------------------------------------- void FundamentalSectorGeometry::precomputeAzimuthalCorrection() { - constexpr double k_TwoPi = 2.0 * 3.14159265358979323846; + constexpr double k_Pi = 3.14159265358979323846; + constexpr double k_TwoPi = 2.0 * k_Pi; + + if(m_Vertices.size() < 2 || m_BoundaryNormals.empty()) + { + // No correction possible -- use identity mapping + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + m_AzimuthalCorrectionTable[i] = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + } + return; + } + + // Step 1: Sample the boundary distance d(rho) at each azimuthal angle. + // For each sampled angle, rotate a reference direction around the barycenter + // by that angle and compute the radial distance to the boundary. + // + // Instead of doing full polarCoordinates (expensive), we directly compute + // the boundary distance: for each angle, create a direction at a small + // offset from the barycenter, then measure how far the boundary is. + + // Reference direction in the tangent plane at the barycenter + Vec3 ref = {0.0, 0.0, 1.0}; + if(std::abs(vecDot(ref, m_Barycenter)) > 0.99) + { + ref = {1.0, 0.0, 0.0}; + } + double refDotCenter = vecDot(ref, m_Barycenter); + Vec3 rx = vecNormalize({ref[0] - refDotCenter * m_Barycenter[0], + ref[1] - refDotCenter * m_Barycenter[1], + ref[2] - refDotCenter * m_Barycenter[2]}); + Vec3 ry = vecNormalize(vecCross(m_Barycenter, rx)); + + // For each sampled angle, compute the angular distance from barycenter to boundary + std::array boundaryDist = {}; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) { - m_AzimuthalCorrectionTable[i] = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + double angle = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + double cosA = std::cos(angle); + double sinA = std::sin(angle); + + // Direction in the tangent plane at this azimuth + Vec3 tangentDir = {cosA * rx[0] + sinA * ry[0], + cosA * rx[1] + sinA * ry[1], + cosA * rx[2] + sinA * ry[2]}; + + // Create a test direction slightly away from barycenter in this tangent direction + // We use a small angle offset (e.g., 0.01 radians) to stay in the linear regime + constexpr double k_SmallAngle = 0.01; + Vec3 testDir = vecNormalize({m_Barycenter[0] + k_SmallAngle * tangentDir[0], + m_Barycenter[1] + k_SmallAngle * tangentDir[1], + m_Barycenter[2] + k_SmallAngle * tangentDir[2]}); + + // Compute the boundary distance at this azimuth using the same algorithm as polarCoordinates + Vec3 gcNormal = vecNormalize(vecCross(m_Barycenter, testDir)); + double distMax = k_Pi; // default large distance + + // Handle degenerate gcNormal (testDir ~= barycenter) + double gcLen = std::sqrt(gcNormal[0] * gcNormal[0] + gcNormal[1] * gcNormal[1] + gcNormal[2] * gcNormal[2]); + if(gcLen < 1.0e-10) + { + boundaryDist[i] = 1.0; + continue; + } + + for(const auto& normal : m_BoundaryNormals) + { + Vec3 bp = vecNormalize(vecCross(normal, gcNormal)); + // Choose the intersection on the side of the barycenter + if(vecDot(testDir, bp) < 0.0) + { + bp = vecNeg(bp); + } + double d = vecAngle(m_Barycenter, bp); + if(d > 1.0e-10) + { + distMax = std::min(distMax, d); + } + } + boundaryDist[i] = distMax; + } + + // Step 2: Compute vertex azimuths and assign equal hue sectors + size_t nVerts = m_Vertices.size(); + + // Compute the azimuthal angle of each vertex relative to the barycenter + std::vector vertexAngles(nVerts); + for(size_t v = 0; v < nVerts; v++) + { + double hDotCenter = vecDot(m_Vertices[v], m_Barycenter); + Vec3 dv = {m_Vertices[v][0] - hDotCenter * m_Barycenter[0], + m_Vertices[v][1] - hDotCenter * m_Barycenter[1], + m_Vertices[v][2] - hDotCenter * m_Barycenter[2]}; + vertexAngles[v] = std::fmod(std::atan2(vecDot(ry, dv), vecDot(rx, dv)) + k_TwoPi, k_TwoPi); + } + + // Sort vertex angles + std::vector sortIdx(nVerts); + std::iota(sortIdx.begin(), sortIdx.end(), 0); + std::sort(sortIdx.begin(), sortIdx.end(), [&](size_t a, size_t b) { return vertexAngles[a] < vertexAngles[b]; }); + std::vector sortedAngles(nVerts); + for(size_t i = 0; i < nVerts; i++) + { + sortedAngles[i] = vertexAngles[sortIdx[i]]; + } + + // Step 3: Build the weighted CDF with boundary distance weighting + // Weight each angular sample by d(rho) -- this is the core of the paper's + // hue speed function. Directions where the boundary is farther get more hue space. + std::array weights = {}; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + weights[i] = boundaryDist[i]; // weight by boundary distance + } + + // Normalize within each vertex sector so each sector gets exactly (2*pi / nVerts) + double sectorSize = k_TwoPi / static_cast(nVerts); + for(size_t s = 0; s < nVerts; s++) + { + double sectorStart = sortedAngles[s]; + double sectorEnd = (s + 1 < nVerts) ? sortedAngles[s + 1] : sortedAngles[0] + k_TwoPi; + + // Find indices in this sector + double sectorSum = 0.0; + size_t count = 0; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + double angle = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + // Check if angle is in this sector (handle wrap-around) + bool inSector = false; + if(sectorEnd <= k_TwoPi) + { + inSector = (angle >= sectorStart && angle < sectorEnd); + } + else + { + inSector = (angle >= sectorStart || angle < std::fmod(sectorEnd, k_TwoPi)); + } + if(inSector) + { + sectorSum += weights[i]; + count++; + } + } + + // Normalize this sector's weights so they sum to sectorSize + if(sectorSum > 1.0e-10 && count > 0) + { + double scale = sectorSize / sectorSum; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + double angle = static_cast(i) / static_cast(k_AzimuthalTableSize) * k_TwoPi; + bool inSector = false; + if(sectorEnd <= k_TwoPi) + { + inSector = (angle >= sectorStart && angle < sectorEnd); + } + else + { + inSector = (angle >= sectorStart || angle < std::fmod(sectorEnd, k_TwoPi)); + } + if(inSector) + { + weights[i] *= scale; + } + } + } + } + + // Step 4: Cumulative sum -> correction table + // The CDF maps raw angle to corrected angle + double cumSum = 0.0; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + cumSum += weights[i]; + m_AzimuthalCorrectionTable[i] = cumSum; + } + + // Normalize so the total is exactly 2*pi + if(cumSum > 1.0e-10) + { + double scale = k_TwoPi / cumSum; + for(size_t i = 0; i < k_AzimuthalTableSize; i++) + { + m_AzimuthalCorrectionTable[i] *= scale; + } } } diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp index 6bd53f4..44779f2 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -197,9 +197,11 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& auto [radius, rho] = m_Sector.polarCoordinates(direction); // 2. Hue from azimuthal angle - // rho is in [0, 2*pi) -- normalize to [0, 1) for HSL conversion - // Then apply Gaussian CDF correction to expand yellow/cyan/magenta regions - double hue = correctHue(rho / k_TwoPi); + // First apply boundary-distance-weighted azimuthal correction to smooth + // the transitions between boundary zones and equalize vertex hue sectors. + // Then apply Gaussian CDF correction to expand yellow/cyan/magenta regions. + double rhoCorrected = m_Sector.correctAzimuthalAngle(rho); + double hue = correctHue(rhoCorrected / k_TwoPi); // 3. Lightness from radial distance via gray gradient blending // @@ -266,7 +268,8 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& if(inSupergroup) { auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); - hue = correctHue(sgRho / k_TwoPi); + double sgRhoCorrected = m_SupergroupSector->correctAzimuthalAngle(sgRho); + hue = correctHue(sgRhoCorrected / k_TwoPi); double rEff = std::pow(sgRadius, 0.35); double r = 1.0 - rEff / 2.0; computeColorFromSphere(r, k_GrayValueWhite); diff --git a/Source/Test/FundamentalSectorGeometryTest.cpp b/Source/Test/FundamentalSectorGeometryTest.cpp index 15dde79..0e5ff3c 100644 --- a/Source/Test/FundamentalSectorGeometryTest.cpp +++ b/Source/Test/FundamentalSectorGeometryTest.cpp @@ -161,12 +161,16 @@ TEST_CASE("ebsdlib::FundamentalSectorGeometry::CorrectAzimuthalAngle", "[EbsdLib { auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); - SECTION("Identity mapping: output approximates input") + SECTION("Correction produces monotonically increasing output") { - // With the identity precomputation, corrected angle should be close to input - double rho = 1.5; - double corrected = sector.correctAzimuthalAngle(rho); - REQUIRE(corrected == Approx(rho).margin(0.02)); + // The corrected angle should increase monotonically with input angle + double prev = 0.0; + for(double rho = 0.01; rho < 2.0 * M_PI - 0.01; rho += 0.05) + { + double corrected = sector.correctAzimuthalAngle(rho); + REQUIRE(corrected >= prev - 0.01); // monotonic (with small tolerance) + prev = corrected; + } } SECTION("Result is in [0, 2*pi)") From ea7451ac39c1e818d588c47ded560484f2194b2f Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Tue, 24 Mar 2026 12:16:14 -0400 Subject: [PATCH 14/16] fix: use orix/MTEX polar coordinate convention to eliminate gray bands The sharp color transition lines ("gray bands") radiating from the center were caused by the radius computation producing non-monotonic values along radial paths. When the "closest boundary" switched between different boundary normals, the radius would dip, creating visible seams. Fixed by adopting the orix/MTEX polar coordinate algorithm: - Uses angle(-v, bp) / angle(-center, bp) instead of angle(center, h) / angle(center, bp) - Cross product order: v.cross(center) then gcN.cross(normal) - Convention: radius=1 at center (white), radius=0 at boundary (saturated) - This produces monotonically varying radius along any radial path, eliminating the gray band artifacts Updated NolzeHielscherColorKey to use the new convention: - White center: r = 0.5 + radius/2 maps [1,0] -> [1.0, 0.5] - No power curve needed since the orix formula already distributes colors well Co-Authored-By: Claude Opus 4.6 (1M context) --- .../Utilities/FundamentalSectorGeometry.cpp | 54 ++++++++++--------- .../Utilities/NolzeHielscherColorKey.cpp | 14 ++--- Source/Test/FundamentalSectorGeometryTest.cpp | 21 ++++---- 3 files changed, 47 insertions(+), 42 deletions(-) diff --git a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp index 9b9c708..3e5d9ef 100644 --- a/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp +++ b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp @@ -131,52 +131,56 @@ std::pair FundamentalSectorGeometry::polarCoordinates(const Vec3 constexpr double k_Pi = 3.14159265358979323846; // Singularity guard: if h is at the barycenter + // Convention: radius=1 at center, 0 at boundary double angleToCenter = vecAngle(h, m_Barycenter); if(angleToCenter < 1.0e-10) { - return {0.0, 0.0}; + return {1.0, 0.0}; } // ------------------------------------------------------------------- // RADIUS: normalized distance from center to boundary // ------------------------------------------------------------------- - // Normal to the great circle through center and h - Vec3 gcNormal = vecNormalize(vecCross(m_Barycenter, h)); + // Algorithm (from orix/MTEX polarCoordinates): + // The great circle plane containing both h and center has normal + // gcN = normalize(v.cross(center)). + // For each boundary normal N_j, the intersection is: + // bp = normalize(gcN.cross(N_j)) + // The radius is: min over j of angle(-v, bp) / angle(-center, bp) + // + // This gives radius=0 at center, radius=1 at boundary, + // and increases monotonically along any radial direction. + // ------------------------------------------------------------------- + Vec3 hNeg = vecNeg(h); + Vec3 centerNeg = vecNeg(m_Barycenter); + + // Normal to the great circle through h and center + // NOTE: order is h cross center (same as orix: v.cross(center)) + Vec3 gcNormal = vecNormalize(vecCross(h, m_Barycenter)); double radius = std::numeric_limits::infinity(); - double distCenterH = angleToCenter; // angle from center to h for(const auto& normal : m_BoundaryNormals) { - // Intersection of the great circle (center,h) with boundary plane N_j. - // The two great circles intersect at: P = normalize(cross(N_j, gcNormal)) and -P. - Vec3 bp = vecNormalize(vecCross(normal, gcNormal)); - - // Choose the intersection point on the same side as h relative to center. - // We want the point P such that the arc center->P passes through or near h. - // Test: if dot(h, bp) < 0, flip to -bp (choose the hemisphere containing h). - if(vecDot(h, bp) < 0.0) - { - bp = vecNeg(bp); - } + // Intersection of great circles: gcN cross N_j + Vec3 bp = vecNormalize(vecCross(gcNormal, normal)); - double distCenterBp = vecAngle(m_Barycenter, bp); + // Compute ratio using antipodal distances + // This naturally selects the correct intersection point + double distNegH = vecAngle(hNeg, bp); + double distNegCenter = vecAngle(centerNeg, bp); double ratio; - if(distCenterBp < 1.0e-10) + if(distNegCenter < 1.0e-10) { - // Center is on this boundary -- boundary is at distance 0, but the sector - // is on the positive side. This shouldn't happen for well-defined sectors. - // Treat as no constraint. - continue; + ratio = 1.0; } else { - ratio = distCenterH / distCenterBp; + ratio = distNegH / distNegCenter; } - // Handle NaN - if(std::isnan(ratio)) + if(std::isnan(ratio) || std::isinf(ratio)) { ratio = 1.0; } @@ -186,7 +190,7 @@ std::pair FundamentalSectorGeometry::polarCoordinates(const Vec3 if(std::isinf(radius)) { - radius = 0.0; + radius = 1.0; } radius = std::clamp(radius, 0.0, 1.0); diff --git a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp index 44779f2..642551c 100644 --- a/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -255,10 +255,9 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& if(m_Sector.colorKeyMode() == "standard" || m_Sector.colorKeyMode() == "impossible") { // Standard: white center only - // Map radius [0,1] -> sphere parameter [1.0, 0.5] - // Use sqrt(radius) to accelerate transition: shrinks white center, expands saturated region - double rEff = std::pow(radius, 0.35); - double r = 1.0 - rEff / 2.0; + // Radius convention: 1 at center, 0 at boundary + // Map to sphere parameter: center(r=1) -> 1.0 (white), boundary(r=0) -> 0.5 (saturated) + double r = 0.5 + radius / 2.0; computeColorFromSphere(r, k_GrayValueWhite); } else if(m_Sector.colorKeyMode() == "extended" && m_SupergroupSector) @@ -270,13 +269,14 @@ NolzeHielscherColorKey::Vec3 NolzeHielscherColorKey::direction2Color(const Vec3& auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); double sgRhoCorrected = m_SupergroupSector->correctAzimuthalAngle(sgRho); hue = correctHue(sgRhoCorrected / k_TwoPi); - double rEff = std::pow(sgRadius, 0.35); - double r = 1.0 - rEff / 2.0; + // White center half: center(r=1)->1.0(white), boundary(r=0)->0.5(saturated) + double r = 0.5 + sgRadius / 2.0; computeColorFromSphere(r, k_GrayValueWhite); } else { - double rEff = std::pow(radius, 0.35); + // Black center half: center(r=1)->0.0(black), boundary(r=0)->0.5(saturated) + double rEff = radius; double r = rEff / 2.0; computeColorFromSphere(r, k_GrayValueBlack); } diff --git a/Source/Test/FundamentalSectorGeometryTest.cpp b/Source/Test/FundamentalSectorGeometryTest.cpp index 0e5ff3c..91c6855 100644 --- a/Source/Test/FundamentalSectorGeometryTest.cpp +++ b/Source/Test/FundamentalSectorGeometryTest.cpp @@ -45,32 +45,33 @@ TEST_CASE("ebsdlib::FundamentalSectorGeometry::PolarCoordinates", "[EbsdLib][Fun { auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); - SECTION("At barycenter: radius = 0") + SECTION("At barycenter: radius = 1 (center of sector)") { auto center = sector.barycenter(); auto [radius, rho] = sector.polarCoordinates(center); - REQUIRE(radius == Approx(0.0).margin(1e-4)); + // Convention: radius=1 at center, 0 at boundary (orix/MTEX convention) + REQUIRE(radius == Approx(1.0).margin(1e-4)); } - SECTION("At vertex [001]: radius near 1") + SECTION("At vertex [001]: radius near 0 (on boundary)") { Vec3 v001 = {0.0, 0.0, 1.0}; auto [radius, rho] = sector.polarCoordinates(v001); - REQUIRE(radius == Approx(1.0).margin(0.05)); + REQUIRE(radius == Approx(0.0).margin(0.05)); } - SECTION("At vertex [101]: radius near 1") + SECTION("At vertex [101]: radius near 0 (on boundary)") { Vec3 v101 = normalize({1.0, 0.0, 1.0}); auto [radius, rho] = sector.polarCoordinates(v101); - REQUIRE(radius == Approx(1.0).margin(0.05)); + REQUIRE(radius == Approx(0.0).margin(0.05)); } - SECTION("At vertex [111]: radius near 1") + SECTION("At vertex [111]: radius near 0 (on boundary)") { Vec3 v111 = normalize({1.0, 1.0, 1.0}); auto [radius, rho] = sector.polarCoordinates(v111); - REQUIRE(radius == Approx(1.0).margin(0.05)); + REQUIRE(radius == Approx(0.0).margin(0.05)); } SECTION("Radius is in [0, 1] for interior point") @@ -95,11 +96,11 @@ TEST_CASE("ebsdlib::FundamentalSectorGeometry::EdgeCases", "[EbsdLib][Fundamenta { auto sector = ebsdlib::FundamentalSectorGeometry::cubicHigh(); - SECTION("Direction exactly at barycenter returns radius = 0") + SECTION("Direction exactly at barycenter returns radius = 1 (center)") { auto center = sector.barycenter(); auto [radius, rho] = sector.polarCoordinates(center); - REQUIRE(radius == Approx(0.0).margin(1e-6)); + REQUIRE(radius == Approx(1.0).margin(1e-6)); } SECTION("isInside returns true for interior, false for exterior") From 23e606c429a7bf26524d2495bf6f6b3846bf6450 Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Fri, 3 Apr 2026 13:35:11 -0400 Subject: [PATCH 15/16] feat: add GriddedColorKey decorator for MTEX-style legend rendering Implements the Decorator pattern: GriddedColorKey wraps any IColorKey with grid-based flat shading. Colors are precomputed at regular angular intervals (default 1 degree, matching MTEX) and looked up via nearest- grid-point snapping. This produces flat-colored patches that hide C1 discontinuities in the underlying color function. Adds LegendRenderMode enum and setLegendRenderMode() convenience method to LaueOps for switching between PerPixel (EbsdLib default) and GridInterpolated (MTEX-style) rendering. Co-Authored-By: Claude Opus 4.6 (1M context) --- Source/Apps/generate_ipf_legends.cpp | 20 ++- Source/EbsdLib/LaueOps/LaueOps.cpp | 27 +++ Source/EbsdLib/LaueOps/LaueOps.h | 10 ++ Source/EbsdLib/Utilities/GriddedColorKey.cpp | 106 ++++++++++++ Source/EbsdLib/Utilities/GriddedColorKey.hpp | 79 +++++++++ Source/EbsdLib/Utilities/SourceList.cmake | 2 + Source/Test/CMakeLists.txt | 1 + Source/Test/GriddedColorKeyTest.cpp | 172 +++++++++++++++++++ 8 files changed, 416 insertions(+), 1 deletion(-) create mode 100644 Source/EbsdLib/Utilities/GriddedColorKey.cpp create mode 100644 Source/EbsdLib/Utilities/GriddedColorKey.hpp create mode 100644 Source/Test/GriddedColorKeyTest.cpp diff --git a/Source/Apps/generate_ipf_legends.cpp b/Source/Apps/generate_ipf_legends.cpp index c344f6c..9c5d768 100644 --- a/Source/Apps/generate_ipf_legends.cpp +++ b/Source/Apps/generate_ipf_legends.cpp @@ -18,6 +18,7 @@ #include "EbsdLib/Utilities/ColorTable.h" #include "EbsdLib/Utilities/EbsdStringUtils.hpp" #include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/GriddedColorKey.hpp" #include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" #include "EbsdLib/Utilities/TSLColorKey.hpp" #include "EbsdLib/Utilities/TiffWriter.h" @@ -374,6 +375,23 @@ void GenerateNolzeHielscherLegends(int imageDim) result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); std::cout << ops.getSymmetryName() << " NH Triangle Result: " << result.first << ": " << result.second << std::endl; + // Set to grid-interpolated mode (MTEX-style rendering) + ops.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 1.0); + + // Generate gridded full-circle legend + legend = ops.generateIPFTriangleLegend(imageDim, true); + ss.str(""); + ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_GRIDDED_FULL.tiff"; + result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + std::cout << ops.getSymmetryName() << " NH Gridded Full Result: " << result.first << ": " << result.second << std::endl; + + // Generate gridded triangle-only legend + legend = ops.generateIPFTriangleLegend(imageDim, false); + ss.str(""); + ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_GRIDDED.tiff"; + result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + std::cout << ops.getSymmetryName() << " NH Gridded Triangle Result: " << result.first << ": " << result.second << std::endl; + // Reset to TSL for subsequent operations ops.setColorKey(std::make_shared()); } @@ -393,7 +411,7 @@ int main(int argc, char* argv[]) } std::stringstream ss; - int imageDim = 512; + int imageDim = 1500; { TrigonalOps ops; auto legend = ops.generateIPFTriangleLegend(imageDim, true); diff --git a/Source/EbsdLib/LaueOps/LaueOps.cpp b/Source/EbsdLib/LaueOps/LaueOps.cpp index 9b82b21..7c8b436 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.cpp +++ b/Source/EbsdLib/LaueOps/LaueOps.cpp @@ -51,6 +51,7 @@ #include "EbsdLib/Orientation/Quaternion.hpp" #include "EbsdLib/Utilities/ColorTable.h" #include "EbsdLib/Utilities/ComputeStereographicProjection.h" +#include "EbsdLib/Utilities/GriddedColorKey.hpp" #include "EbsdLib/Utilities/TSLColorKey.hpp" #include // for std::max @@ -113,6 +114,32 @@ ebsdlib::IColorKey::Pointer LaueOps::getColorKey() const return m_ColorKey; } +// ----------------------------------------------------------------------------- +void LaueOps::setLegendRenderMode(ebsdlib::LegendRenderMode mode, double gridResolutionDeg) +{ + if(mode == ebsdlib::LegendRenderMode::GridInterpolated) + { + // Wrap the current color key with a GriddedColorKey if not already wrapped + auto currentKey = m_ColorKey; + // If already gridded, unwrap first to avoid double-wrapping + auto griddedKey = std::dynamic_pointer_cast(currentKey); + if(griddedKey) + { + currentKey = griddedKey->innerKey(); + } + m_ColorKey = std::make_shared(currentKey, gridResolutionDeg); + } + else + { + // PerPixel mode: unwrap if currently gridded + auto griddedKey = std::dynamic_pointer_cast(m_ColorKey); + if(griddedKey) + { + m_ColorKey = griddedKey->innerKey(); + } + } +} + // ----------------------------------------------------------------------------- std::string LaueOps::FZTypeToString(const FZType value) { diff --git a/Source/EbsdLib/LaueOps/LaueOps.h b/Source/EbsdLib/LaueOps/LaueOps.h index e0084c6..a71b116 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.h +++ b/Source/EbsdLib/LaueOps/LaueOps.h @@ -46,6 +46,7 @@ #include "EbsdLib/Orientation/OrientationFwd.hpp" #include "EbsdLib/Orientation/Quaternion.hpp" #include "EbsdLib/Orientation/Rodrigues.hpp" +#include "EbsdLib/Utilities/GriddedColorKey.hpp" #include "EbsdLib/Utilities/IColorKey.hpp" #include "EbsdLib/Utilities/PoleFigureUtilities.h" #include "EbsdLib/Utilities/TSLColorKey.hpp" @@ -302,6 +303,15 @@ class EbsdLib_EXPORT LaueOps */ ebsdlib::IColorKey::Pointer getColorKey() const; + /** + * @brief Set the legend rendering mode. + * PerPixel: exact color at every pixel (default, current behavior) + * GridInterpolated: MTEX-style flat-shaded grid cells at the given resolution + * @param mode The rendering mode to use + * @param gridResolutionDeg Grid cell size in degrees (only used for GridInterpolated mode) + */ + void setLegendRenderMode(ebsdlib::LegendRenderMode mode, double gridResolutionDeg = 1.0); + /** * @brief generateRodriguesColor Generates an RGB Color from a Rodrigues Vector * @param r1 First component of the Rodrigues Vector diff --git a/Source/EbsdLib/Utilities/GriddedColorKey.cpp b/Source/EbsdLib/Utilities/GriddedColorKey.cpp new file mode 100644 index 0000000..7354a14 --- /dev/null +++ b/Source/EbsdLib/Utilities/GriddedColorKey.cpp @@ -0,0 +1,106 @@ +#include "EbsdLib/Utilities/GriddedColorKey.hpp" + +#include +#include + +namespace ebsdlib +{ + +namespace +{ +constexpr double k_Pi = 3.14159265358979323846; +constexpr double k_HalfPi = k_Pi / 2.0; +constexpr double k_DegToRad = k_Pi / 180.0; +} // namespace + +GriddedColorKey::GriddedColorKey(IColorKey::Pointer innerKey, double resolutionDeg) +: m_InnerKey(std::move(innerKey)) +, m_ResolutionDeg(resolutionDeg) +, m_ResolutionRad(resolutionDeg * k_DegToRad) +{ + // Grid covers eta in [0, pi] (180 degrees) and chi in [0, pi/2] (90 degrees) + // This covers all possible Laue group SSTs + m_EtaSteps = static_cast(std::ceil(180.0 / resolutionDeg)) + 1; + m_ChiSteps = static_cast(std::ceil(90.0 / resolutionDeg)) + 1; + precomputeGrid(); +} + +void GriddedColorKey::precomputeGrid() +{ + m_Grid.resize(m_EtaSteps); + + for(int ei = 0; ei < m_EtaSteps; ei++) + { + m_Grid[ei].resize(m_ChiSteps); + double eta = static_cast(ei) * m_ResolutionRad; + + for(int ci = 0; ci < m_ChiSteps; ci++) + { + double chi = static_cast(ci) * m_ResolutionRad; + + // Convert spherical to Cartesian direction and compute color + // via the inner key's Vec3 overload + double sinChi = std::sin(chi); + double cosChi = std::cos(chi); + Vec3 dir = {sinChi * std::cos(eta), sinChi * std::sin(eta), cosChi}; + + m_Grid[ei][ci] = m_InnerKey->direction2Color(dir); + } + } +} + +GriddedColorKey::Vec3 GriddedColorKey::lookupGrid(double eta, double chi) const +{ + // Map to grid indices via nearest-neighbor snapping + int ei = static_cast(std::round(eta / m_ResolutionRad)); + int ci = static_cast(std::round(chi / m_ResolutionRad)); + + // Clamp to grid bounds + ei = std::clamp(ei, 0, m_EtaSteps - 1); + ci = std::clamp(ci, 0, m_ChiSteps - 1); + + return m_Grid[ei][ci]; +} + +GriddedColorKey::Vec3 GriddedColorKey::direction2Color(const Vec3& direction) const +{ + // Convert direction to (eta, chi) and look up from grid + double chi = std::acos(std::clamp(direction[2], -1.0, 1.0)); + double eta = std::atan2(direction[1], direction[0]); + if(eta < 0.0) + { + eta += 2.0 * k_Pi; + } + return lookupGrid(eta, chi); +} + +GriddedColorKey::Vec3 GriddedColorKey::direction2Color(double eta, double chi, const Vec3& angleLimits) const +{ + // Snap eta and chi to nearest grid point (flat shading). + // Instead of computing the exact color at (eta, chi), + // we return the precomputed color at the nearest grid point. + // This produces flat-colored patches like MTEX's surf() rendering. + double etaPositive = eta; + if(etaPositive < 0.0) + { + etaPositive += 2.0 * k_Pi; + } + return lookupGrid(etaPositive, chi); +} + +std::string GriddedColorKey::name() const +{ + return m_InnerKey->name() + " (gridded)"; +} + +IColorKey::Pointer GriddedColorKey::innerKey() const +{ + return m_InnerKey; +} + +double GriddedColorKey::resolutionDeg() const +{ + return m_ResolutionDeg; +} + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/GriddedColorKey.hpp b/Source/EbsdLib/Utilities/GriddedColorKey.hpp new file mode 100644 index 0000000..7ae75e4 --- /dev/null +++ b/Source/EbsdLib/Utilities/GriddedColorKey.hpp @@ -0,0 +1,79 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" +#include "EbsdLib/Utilities/IColorKey.hpp" + +#include +#include +#include + +namespace ebsdlib +{ + +/** + * @brief Rendering mode for IPF legend generation. + */ +enum class LegendRenderMode +{ + PerPixel, ///< Compute exact color at every pixel (EbsdLib default) + GridInterpolated ///< Sample at grid points, flat-shade cells (MTEX-style) +}; + +/** + * @brief Decorator that wraps any IColorKey with grid-based flat shading. + * + * On construction, precomputes colors at a regular grid of (eta, chi) sample + * points by calling the inner color key. When direction2Color() is called, + * it snaps the direction to the nearest grid point and returns the precomputed + * color (flat shading), producing smooth color patches that hide C1 + * discontinuities in the underlying color function. + * + * This replicates the MTEX rendering approach where colors are sampled at + * ~1-degree intervals and rendered as flat-colored quadrilateral patches. + * + * Usage: + * auto nhKey = std::make_shared(sector); + * auto gridKey = std::make_shared(nhKey, 1.0); + * ops.setColorKey(gridKey); + */ +class EbsdLib_EXPORT GriddedColorKey : public IColorKey +{ +public: + /** + * @brief Construct a grid-decorated color key. + * @param innerKey The underlying color key to sample from + * @param resolutionDeg Grid cell size in degrees (default 1.0, matching MTEX) + */ + explicit GriddedColorKey(IColorKey::Pointer innerKey, double resolutionDeg = 1.0); + ~GriddedColorKey() override = default; + + Vec3 direction2Color(const Vec3& direction) const override; + Vec3 direction2Color(double eta, double chi, const Vec3& angleLimits) const override; + std::string name() const override; + + /** + * @brief Get the underlying (unwrapped) color key. + */ + IColorKey::Pointer innerKey() const; + + /** + * @brief Get the grid resolution in degrees. + */ + double resolutionDeg() const; + +private: + IColorKey::Pointer m_InnerKey; + double m_ResolutionRad; // grid cell size in radians + + // Precomputed color grid: m_Grid[etaIdx][chiIdx] = RGB color + // Covers eta in [0, pi] and chi in [0, pi/2] to handle all Laue groups + std::vector> m_Grid; + int m_EtaSteps; + int m_ChiSteps; + double m_ResolutionDeg; + + void precomputeGrid(); + Vec3 lookupGrid(double eta, double chi) const; +}; + +} // namespace ebsdlib diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index cffc4ed..74ef2cc 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -28,6 +28,7 @@ set(EbsdLib_${DIR_NAME}_HDRS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.hpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/NolzeHielscherColorKey.hpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/GriddedColorKey.hpp ) set(EbsdLib_${DIR_NAME}_SRCS @@ -47,6 +48,7 @@ set(EbsdLib_${DIR_NAME}_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/FundamentalSectorGeometry.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/TSLColorKey.cpp ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/NolzeHielscherColorKey.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/GriddedColorKey.cpp ) # # QT5_WRAP_CPP( EbsdLib_Generated_MOC_SRCS ${EbsdLib_Utilities_MOC_HDRS} ) # set_source_files_properties( ${EbsdLib_Generated_MOC_SRCS} PROPERTIES HEADER_FILE_ONLY TRUE) diff --git a/Source/Test/CMakeLists.txt b/Source/Test/CMakeLists.txt index 6b6a017..74ffa8e 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -58,6 +58,7 @@ set(EbsdLib_UnitTest_SRCS ${EbsdLibProj_SOURCE_DIR}/Source/Test/FundamentalSectorGeometryTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/TSLColorKeyTest.cpp ${EbsdLibProj_SOURCE_DIR}/Source/Test/NolzeHielscherColorKeyTest.cpp + ${EbsdLibProj_SOURCE_DIR}/Source/Test/GriddedColorKeyTest.cpp ) diff --git a/Source/Test/GriddedColorKeyTest.cpp b/Source/Test/GriddedColorKeyTest.cpp new file mode 100644 index 0000000..3608cf8 --- /dev/null +++ b/Source/Test/GriddedColorKeyTest.cpp @@ -0,0 +1,172 @@ +#include + +#include "EbsdLib/LaueOps/CubicOps.h" +#include "EbsdLib/LaueOps/LaueOps.h" +#include "EbsdLib/Utilities/ColorTable.h" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/GriddedColorKey.hpp" +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" +#include "EbsdLib/Utilities/TSLColorKey.hpp" + +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +TEST_CASE("ebsdlib::GriddedColorKey::BasicProperties", "[EbsdLib][GriddedColorKey]") +{ + auto tslKey = std::make_shared(); + auto gridKey = std::make_shared(tslKey, 2.0); + + SECTION("Name includes gridded suffix") + { + REQUIRE(gridKey->name() == "TSL (gridded)"); + } + + SECTION("Inner key is accessible") + { + REQUIRE(gridKey->innerKey()->name() == "TSL"); + } + + SECTION("Resolution is stored correctly") + { + REQUIRE(gridKey->resolutionDeg() == Approx(2.0)); + } +} + +TEST_CASE("ebsdlib::GriddedColorKey::FlatShading", "[EbsdLib][GriddedColorKey]") +{ + auto nhKey = std::make_shared(ebsdlib::FundamentalSectorGeometry::cubicHigh()); + auto gridKey = std::make_shared(nhKey, 2.0); // coarse 2-degree grid + + SECTION("Nearby points within same grid cell produce identical colors") + { + // Two points that are less than 2 degrees apart should snap to the same grid cell + double eta1 = 0.2; + double chi1 = 0.3; + double eta2 = 0.2 + 0.01; // ~0.6 degrees apart + double chi2 = 0.3 + 0.01; + + std::array limits = {0.0, M_PI / 4.0, 0.6}; + auto c1 = gridKey->direction2Color(eta1, chi1, limits); + auto c2 = gridKey->direction2Color(eta2, chi2, limits); + + // Should be exactly equal (same grid cell) + REQUIRE(c1[0] == Approx(c2[0]).margin(1e-10)); + REQUIRE(c1[1] == Approx(c2[1]).margin(1e-10)); + REQUIRE(c1[2] == Approx(c2[2]).margin(1e-10)); + } + + SECTION("Points in different grid cells may produce different colors") + { + double eta1 = 0.2; + double eta2 = 0.2 + 0.05; // ~2.9 degrees apart, different cell + + std::array limits = {0.0, M_PI / 4.0, 0.6}; + auto c1 = gridKey->direction2Color(eta1, 0.3, limits); + auto c2 = gridKey->direction2Color(eta2, 0.3, limits); + + // These may or may not differ depending on the color function + // Just verify they are valid + REQUIRE(c1[0] >= 0.0); + REQUIRE(c1[0] <= 1.0); + REQUIRE(c2[0] >= 0.0); + REQUIRE(c2[0] <= 1.0); + } + + SECTION("All outputs are valid RGB") + { + for(double eta = 0.01; eta < M_PI / 4.0 - 0.01; eta += 0.05) + { + double chiMax = std::acos(std::sqrt(1.0 / (2.0 + std::tan(eta) * std::tan(eta)))); + for(double chi = 0.01; chi < chiMax - 0.01; chi += 0.05) + { + std::array limits = {0.0, M_PI / 4.0, chiMax}; + auto [r, g, b] = gridKey->direction2Color(eta, chi, limits); + REQUIRE(r >= 0.0); + REQUIRE(r <= 1.0); + REQUIRE(g >= 0.0); + REQUIRE(g <= 1.0); + REQUIRE(b >= 0.0); + REQUIRE(b <= 1.0); + } + } + } +} + +TEST_CASE("ebsdlib::GriddedColorKey::LaueOpsIntegration", "[EbsdLib][GriddedColorKey]") +{ + auto allOps = ebsdlib::LaueOps::GetAllOrientationOps(); + auto& cubicOps = *allOps[1]; // Cubic_High + + SECTION("Can set gridded color key on LaueOps") + { + auto nhKey = std::make_shared(ebsdlib::FundamentalSectorGeometry::cubicHigh()); + auto gridKey = std::make_shared(nhKey, 1.0); + cubicOps.setColorKey(gridKey); + REQUIRE(cubicOps.getColorKey()->name() == "NolzeHielscher (gridded)"); + + // Generate a legend with the gridded key + auto legend = cubicOps.generateIPFTriangleLegend(64, false); + REQUIRE(legend != nullptr); + REQUIRE(legend->getNumberOfTuples() > 0); + + // Reset + cubicOps.setColorKey(std::make_shared()); + } + + SECTION("IPF colors work with gridded key") + { + auto nhKey = std::make_shared(ebsdlib::FundamentalSectorGeometry::cubicHigh()); + auto gridKey = std::make_shared(nhKey, 1.0); + cubicOps.setColorKey(gridKey); + + double eulers[3] = {0.5, 0.3, 0.2}; + double refDir[3] = {0.0, 0.0, 1.0}; + ebsdlib::Rgb color = cubicOps.generateIPFColor(eulers, refDir, false); + int r = ebsdlib::RgbColor::dRed(color); + int g = ebsdlib::RgbColor::dGreen(color); + int b = ebsdlib::RgbColor::dBlue(color); + REQUIRE((r + g + b) > 0); + + cubicOps.setColorKey(std::make_shared()); + } +} + +TEST_CASE("ebsdlib::GriddedColorKey::SetLegendRenderMode", "[EbsdLib][GriddedColorKey]") +{ + auto allOps = ebsdlib::LaueOps::GetAllOrientationOps(); + auto& cubicOps = *allOps[1]; // Cubic_High + + SECTION("Switch to GridInterpolated mode wraps the color key") + { + cubicOps.setColorKey(std::make_shared()); + cubicOps.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 2.0); + REQUIRE(cubicOps.getColorKey()->name() == "TSL (gridded)"); + } + + SECTION("Switch back to PerPixel mode unwraps the color key") + { + cubicOps.setColorKey(std::make_shared()); + cubicOps.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 2.0); + cubicOps.setLegendRenderMode(ebsdlib::LegendRenderMode::PerPixel); + REQUIRE(cubicOps.getColorKey()->name() == "TSL"); + } + + SECTION("Double-wrapping is prevented") + { + cubicOps.setColorKey(std::make_shared()); + cubicOps.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 2.0); + cubicOps.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 1.0); + // Should still have only one layer of wrapping + REQUIRE(cubicOps.getColorKey()->name() == "TSL (gridded)"); + auto griddedKey = std::dynamic_pointer_cast(cubicOps.getColorKey()); + REQUIRE(griddedKey != nullptr); + REQUIRE(griddedKey->resolutionDeg() == Approx(1.0)); + REQUIRE(griddedKey->innerKey()->name() == "TSL"); + } + + // Reset to default + cubicOps.setColorKey(std::make_shared()); +} From 3874e01d2d97e57e8a9ecf3747a6a8dc3da6031e Mon Sep 17 00:00:00 2001 From: Michael Jackson Date: Fri, 3 Apr 2026 13:41:22 -0400 Subject: [PATCH 16/16] chore: generate gridded NH legends at 2000x2000 with 0.5-degree grid Co-Authored-By: Claude Opus 4.6 (1M context) --- Source/Apps/generate_ipf_legends.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/Source/Apps/generate_ipf_legends.cpp b/Source/Apps/generate_ipf_legends.cpp index 9c5d768..d4d6a7d 100644 --- a/Source/Apps/generate_ipf_legends.cpp +++ b/Source/Apps/generate_ipf_legends.cpp @@ -375,21 +375,22 @@ void GenerateNolzeHielscherLegends(int imageDim) result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); std::cout << ops.getSymmetryName() << " NH Triangle Result: " << result.first << ": " << result.second << std::endl; - // Set to grid-interpolated mode (MTEX-style rendering) - ops.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 1.0); + // Set to grid-interpolated mode (MTEX-style rendering, 0.5 degree grid) + ops.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 0.5); - // Generate gridded full-circle legend - legend = ops.generateIPFTriangleLegend(imageDim, true); + // Generate gridded legends at higher resolution (2000x2000) + constexpr int k_GriddedImageDim = 2000; + legend = ops.generateIPFTriangleLegend(k_GriddedImageDim, true); ss.str(""); ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_GRIDDED_FULL.tiff"; - result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + result = TiffWriter::WriteColorImage(ss.str(), k_GriddedImageDim, k_GriddedImageDim, 3, legend->getPointer(0)); std::cout << ops.getSymmetryName() << " NH Gridded Full Result: " << result.first << ": " << result.second << std::endl; // Generate gridded triangle-only legend - legend = ops.generateIPFTriangleLegend(imageDim, false); + legend = ops.generateIPFTriangleLegend(k_GriddedImageDim, false); ss.str(""); ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_GRIDDED.tiff"; - result = TiffWriter::WriteColorImage(ss.str(), imageDim, imageDim, 3, legend->getPointer(0)); + result = TiffWriter::WriteColorImage(ss.str(), k_GriddedImageDim, k_GriddedImageDim, 3, legend->getPointer(0)); std::cout << ops.getSymmetryName() << " NH Gridded Triangle Result: " << result.first << ": " << result.second << std::endl; // Reset to TSL for subsequent operations