diff --git a/Source/Apps/generate_ipf_legends.cpp b/Source/Apps/generate_ipf_legends.cpp index 705f691..d4d6a7d 100644 --- a/Source/Apps/generate_ipf_legends.cpp +++ b/Source/Apps/generate_ipf_legends.cpp @@ -17,6 +17,10 @@ #include "EbsdLib/Utilities/CanvasUtilities.hpp" #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" #include "EbsdLib/Apps/EbsdLibFileLocations.h" @@ -25,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -324,6 +329,75 @@ 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; + + // Set to grid-interpolated mode (MTEX-style rendering, 0.5 degree grid) + ops.setLegendRenderMode(ebsdlib::LegendRenderMode::GridInterpolated, 0.5); + + // 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(), 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(k_GriddedImageDim, false); + ss.str(""); + ss << k_Output_Dir << "/" << symName << "/" << symName << "_NH_GRIDDED.tiff"; + 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 + ops.setColorKey(std::make_shared()); + } +} + // ----------------------------------------------------------------------------- int main(int argc, char* argv[]) { @@ -338,7 +412,7 @@ int main(int argc, char* argv[]) } std::stringstream ss; - int imageDim = 512; + int imageDim = 1500; { TrigonalOps ops; auto legend = ops.generateIPFTriangleLegend(imageDim, true); @@ -697,5 +771,7 @@ int main(int argc, char* argv[]) GeneratePoleFigures(ops, 1); } + GenerateNolzeHielscherLegends(imageDim); + return 0; } diff --git a/Source/EbsdLib/LaueOps/LaueOps.cpp b/Source/EbsdLib/LaueOps/LaueOps.cpp index 5909bfc..7c8b436 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.cpp +++ b/Source/EbsdLib/LaueOps/LaueOps.cpp @@ -51,6 +51,8 @@ #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 #include @@ -92,11 +94,52 @@ 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; +} + +// ----------------------------------------------------------------------------- +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) { @@ -201,6 +244,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..a71b116 100644 --- a/Source/EbsdLib/LaueOps/LaueOps.h +++ b/Source/EbsdLib/LaueOps/LaueOps.h @@ -46,7 +46,10 @@ #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" namespace ebsdlib { @@ -288,6 +291,27 @@ 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 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 @@ -507,6 +531,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/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/FundamentalSectorGeometry.cpp b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp new file mode 100644 index 0000000..3e5d9ef --- /dev/null +++ b/Source/EbsdLib/Utilities/FundamentalSectorGeometry.cpp @@ -0,0 +1,855 @@ +#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 + // Convention: radius=1 at center, 0 at boundary + double angleToCenter = vecAngle(h, m_Barycenter); + if(angleToCenter < 1.0e-10) + { + return {1.0, 0.0}; + } + + // ------------------------------------------------------------------- + // RADIUS: normalized distance from center to boundary + // ------------------------------------------------------------------- + // 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(); + + for(const auto& normal : m_BoundaryNormals) + { + // Intersection of great circles: gcN cross N_j + Vec3 bp = vecNormalize(vecCross(gcNormal, normal)); + + // 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(distNegCenter < 1.0e-10) + { + ratio = 1.0; + } + else + { + ratio = distNegH / distNegCenter; + } + + if(std::isnan(ratio) || std::isinf(ratio)) + { + ratio = 1.0; + } + + radius = std::min(radius, ratio); + } + + if(std::isinf(radius)) + { + radius = 1.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 +// +// 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_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++) + { + 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; + } + } +} + +// ----------------------------------------------------------------------- +// 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/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/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/NolzeHielscherColorKey.cpp b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp new file mode 100644 index 0000000..642551c --- /dev/null +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.cpp @@ -0,0 +1,298 @@ +#include "EbsdLib/Utilities/NolzeHielscherColorKey.hpp" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.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; +} +/** + * @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 + +// ----------------------------------------------------------------------- +// Constructor +// ----------------------------------------------------------------------- +NolzeHielscherColorKey::NolzeHielscherColorKey(const FundamentalSectorGeometry& sector, double lambdaL, double lambdaS) +: m_Sector(sector) +, 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()); + } + 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; +} + +// ----------------------------------------------------------------------- +// 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; +} + +// ----------------------------------------------------------------------- +// 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) +// ----------------------------------------------------------------------- +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 +// +// 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 +{ + // 1. Get polar coordinates from the fundamental sector geometry + auto [radius, rho] = m_Sector.polarCoordinates(direction); + + // 2. Hue from azimuthal angle + // 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 + // + // 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; + + // 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] + + 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 envelope: peak saturation at th=0.5, reduced at poles + double gray = 1.0 - 2.0 * grayValue * std::abs(th - 0.5); + + // 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 + // 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) + { + bool inSupergroup = m_SupergroupSector->isInside(direction); + + if(inSupergroup) + { + auto [sgRadius, sgRho] = m_SupergroupSector->polarCoordinates(direction); + double sgRhoCorrected = m_SupergroupSector->correctAzimuthalAngle(sgRho); + hue = correctHue(sgRhoCorrected / k_TwoPi); + // 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 + { + // 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); + } + } + + // 5. Convert HSL to RGB + auto rgb = color::hslToRgb(hue, sHsl, 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..3b8fde3 --- /dev/null +++ b/Source/EbsdLib/Utilities/NolzeHielscherColorKey.hpp @@ -0,0 +1,109 @@ +#pragma once + +#include "EbsdLib/EbsdLib.h" +#include "EbsdLib/Utilities/FundamentalSectorGeometry.hpp" +#include "EbsdLib/Utilities/IColorKey.hpp" + +#include +#include +#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); + + /** + * @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 diff --git a/Source/EbsdLib/Utilities/SourceList.cmake b/Source/EbsdLib/Utilities/SourceList.cmake index 0e7fc72..74ef2cc 100644 --- a/Source/EbsdLib/Utilities/SourceList.cmake +++ b/Source/EbsdLib/Utilities/SourceList.cmake @@ -23,6 +23,12 @@ 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 + ${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 + ${EbsdLibProj_SOURCE_DIR}/Source/EbsdLib/${DIR_NAME}/GriddedColorKey.hpp ) set(EbsdLib_${DIR_NAME}_SRCS @@ -39,6 +45,10 @@ 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 + ${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/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 c0754d1..74ffa8e 100644 --- a/Source/Test/CMakeLists.txt +++ b/Source/Test/CMakeLists.txt @@ -54,6 +54,11 @@ 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 + ${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/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)); + } +} diff --git a/Source/Test/FundamentalSectorGeometryTest.cpp b/Source/Test/FundamentalSectorGeometryTest.cpp new file mode 100644 index 0000000..91c6855 --- /dev/null +++ b/Source/Test/FundamentalSectorGeometryTest.cpp @@ -0,0 +1,291 @@ +#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 = 1 (center of sector)") + { + auto center = sector.barycenter(); + auto [radius, rho] = sector.polarCoordinates(center); + // 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 0 (on boundary)") + { + Vec3 v001 = {0.0, 0.0, 1.0}; + auto [radius, rho] = sector.polarCoordinates(v001); + REQUIRE(radius == Approx(0.0).margin(0.05)); + } + + 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(0.0).margin(0.05)); + } + + 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(0.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 = 1 (center)") + { + auto center = sector.barycenter(); + auto [radius, rho] = sector.polarCoordinates(center); + REQUIRE(radius == Approx(1.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("Correction produces monotonically increasing output") + { + // 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)") + { + 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())); + } +} 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()); +} 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()); + } + } +} diff --git a/Source/Test/NolzeHielscherColorKeyTest.cpp b/Source/Test/NolzeHielscherColorKeyTest.cpp new file mode 100644 index 0000000..8d05e54 --- /dev/null +++ b/Source/Test/NolzeHielscherColorKeyTest.cpp @@ -0,0 +1,431 @@ +#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.8); + } + + 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); + } +} + +// --------------------------------------------------------------------------- +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); + } + } +} + +// --------------------------------------------------------------------------- +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.8); + } +} + +// --------------------------------------------------------------------------- +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); + } +} diff --git a/Source/Test/TSLColorKeyTest.cpp b/Source/Test/TSLColorKeyTest.cpp new file mode 100644 index 0000000..441a8ff --- /dev/null +++ b/Source/Test/TSLColorKeyTest.cpp @@ -0,0 +1,234 @@ +#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 +#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); +} + +// --------------------------------------------------------------------------- +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); + } + } +}