Create colors from wavelengths and estimate dominant wavelengths (#504)

Author: facelessuser

coloraide/algebra.py

@@ -1182,6 +1182,61 @@ def pprint(value: float | ArrayLike) -> None:
     print(pretty(value))
 
 
+def ray_line_intersect(
+    a1: VectorLike,
+    a2: VectorLike,
+    b1: VectorLike,
+    b2: VectorLike,
+    abs_tol: float = ATOL
+) -> Vector | None:
+    """Find the intersection of a ray and a line."""
+
+    da = [a - b for a, b in zip(a2, a1)]
+    db = [a - b for a, b in zip(b2, b1)]
+    dp = [a - b for a, b in zip(a1, b1)]
+    dap = [-da[1], da[0]]
+    denom = dot(dap, db, dims=D1)
+    # Perpendicular cases
+    if abs(denom) < abs_tol:  # pragma: no cover
+        return None
+    t = dot(dap, dp, dims=D1) / denom
+    # Intersect
+    i = add(b1, multiply(t, db, dims=SC_D1), dims=D1)
+    # Check if intersection is within bounds
+    if 0 <= t <= 1:
+        return i
+    return None  # pragma: no cover
+
+
+def line_intersect(
+    a1: VectorLike,
+    a2: VectorLike,
+    b1: VectorLike,
+    b2: VectorLike,
+    abs_tol: float = ATOL
+) -> Vector | None:  # pragma: no cover
+    """Find the intersection of of lines."""
+
+    da = [a - b for a, b in zip(a2, a1)]
+    db = [a - b for a, b in zip(b2, b1)]
+    dp = [a - b for a, b in zip(a1, b1)]
+    dap = [-da[1], da[0]]
+    denom = dot(dap, db, dims=D1)
+    # Perpendicular cases
+    if abs(denom) < abs_tol:  # pragma: no cover
+        return None
+    t = dot(dap, dp, dims=D1) / denom
+    # Intersect
+    i = add(b1, multiply(t, db, dims=SC_D1), dims=D1)
+    # Check if intersection is within bounds
+    if (
+        0 <= t <= 1 and
+        0 <= divide(dot(subtract(i, a1, dims=D1), da, dims=D1), dot(da, da, dims=D1), dims=SC) <= 1
+    ):
+        return i
+    return None
+
+
 def all(a: float | ArrayLike) -> bool:  # noqa: A001
     """Return true if all elements are "true"."""
 

coloraide/cmfs.py

@@ -13,10 +13,10 @@ def __init__(self, cmfs: dict[float, Vector], /, interpolator: str = 'cubic', **
         """Initialize and capture attributes of CMF."""
 
         keys = list(cmfs.keys())
-        self.spline = alg.interpolate(list(cmfs.values()), method='sprague')
         self.start = int(keys[0])
         self.end = int(keys[-1])
         self.step = round((self.end - self.start) / (len(keys) - 1))
+        self.spline = alg.interpolate(list(cmfs.values()), method='sprague', domain=[self.start, self.end])
         super().__init__(cmfs, **kwargs)
 
     def __getitem__(self, key: float) -> Vector:  # pragma: no cover
@@ -25,9 +25,7 @@ def __getitem__(self, key: float) -> Vector:  # pragma: no cover
         value = super().get(key)
         if value is not None:
             return value[:]
-
-        factor = (key - self.start) / (self.end - self.start)
-        return self.spline(factor)
+        return self.spline(key)
 
     def xy(self, key: float) -> Vector:  # pragma: no cover
         """Return CMFs data as XY coordinates."""

coloraide/color.py

@@ -19,6 +19,7 @@
 from . import temperature
 from . import util
 from . import algebra as alg
+from . import spectrum
 from .channels import ANGLE_DEG, ANGLE_RAD, ANGLE_GRAD, ANGLE_TURN, ANGLE_NULL
 from .deprecate import warn_deprecated, deprecated
 from itertools import zip_longest as zipl
@@ -1342,6 +1343,51 @@ def contrast(self, color: ColorInput, method: str | None = None) -> float:
         color = self._handle_color_input(color)
         return contrast.contrast(method, self, color)
 
+    def wavelength(
+        self,
+        *,
+        white: VectorLike | None = None,
+        complementary: bool = False
+    ) -> tuple[float, Vector, Vector]:
+        """Get the dominant wavelength."""
+
+        return spectrum.closest_wavelength(
+            self.xy(),
+            white or util.xyz_to_xyY(self.white())[:-1],
+            reverse=complementary
+        )
+
+    @classmethod
+    def from_wavelength(
+        cls,
+        space: str,
+        wavelength: float,
+        *,
+        scale: bool = True,
+        scale_space: str | None = None
+    ) -> Self:
+        """Create a color from a wavelength."""
+
+        xyz = spectrum.wavelength_to_color(wavelength)
+        color = cls('xyz-d65', xyz)
+
+        if scale_space is None:
+            scale_space = 'srgb-linear'
+
+        # Bypass chromatic adaptation to ensure we get the appropriate wavelength into the given white point.
+        # If a scaling is applied, we skip to get to this space, but subsequent conversions will be adapted.
+        if scale and isinstance(cls.CS_MAP[scale_space], RGBish):
+            color._space, color._coords[:-1] = convert.convert(color, scale_space, skip_adaptation=True)
+            color[:-1] = util.rgb_scale(color.coords())
+            # Convert to targeted color space
+            color.convert(space, in_place=True)
+        else:
+            color._space, color._coords[:-1] = convert.convert(color, space, skip_adaptation=True)
+            # Unlikely to get an achromatic, but just in case
+            if color._space.is_polar() and color.is_achromatic():  # pragma: no cover
+                color[color._space.hue_index()] = math.nan  # type: ignore[attr-defined]
+        return color
+
     @overload
     def get(self,
         name: str,

coloraide/convert.py

@@ -119,7 +119,7 @@ def get_convert_chain(
     return chain
 
 
-def convert(color: Color, space: str) -> tuple[Space, Vector]:
+def convert(color: Color, space: str, skip_adaptation: bool = False) -> tuple[Space, Vector]:
     """Convert the color coordinates to the specified space."""
 
     # Grab the convert for the current space to the desired space
@@ -133,15 +133,15 @@ def convert(color: Color, space: str) -> tuple[Space, Vector]:
     # Perform chromatic adaption if needed (a conversion to or from XYZ D65).
     last = color._space
     for a, b, direction, adapt in chain:
-        if direction and adapt:
+        if direction and adapt and not skip_adaptation:
             coords = color.chromatic_adaptation(
                 a.WHITE,
                 b.WHITE,
                 coords
             )
 
         coords = b.from_base(coords) if direction else a.to_base(coords)
-        if not direction and adapt:
+        if not direction and adapt and not skip_adaptation:
             coords = color.chromatic_adaptation(
                 a.WHITE,
                 b.WHITE,

coloraide/spectrum.py

@@ -0,0 +1,157 @@
+"""
+Handle spectral related things.
+
+- Calculations colors from/to wavelengths.
+"""
+from __future__ import annotations
+import math
+from . import algebra as alg
+from . import util
+from .types import Vector, VectorLike
+from . import cat
+from . import cmfs
+
+WHITE = cat.WHITES['2deg']['E']
+LOCUS_START = 360
+LOCUS_END = 780
+LOCUS_STEP = 1
+
+
+def xy_to_angle(xy: VectorLike, white: VectorLike, offset: float = 0.0, invert: bool = False) -> float:
+    """
+    Translate xy to an angle with white being the origin.
+
+    If offset is provided, make the angle relative to the offset angle.
+
+    If invert is requested, we want to rotate the xy point 180 degrees.
+    """
+
+    norm = alg.subtract(xy, white, dims=alg.D1)
+    if invert:
+        norm = alg.multiply(norm, -1, dims=alg.D1_SC)
+    if offset:
+        angle = (alg.rect_to_polar(*norm)[1] - offset) % 360
+        if angle < 1e-12:
+            angle = 360.0
+    else:
+        angle = alg.rect_to_polar(*norm)[1]
+    return angle
+
+
+def closest_wavelength(
+    xy: VectorLike,
+    white: VectorLike = WHITE,
+    reverse: bool = False,
+    closest: bool = True
+) -> tuple[float, Vector, Vector]:
+    """
+    Get the closest dominant wavelength.
+
+    Both intersections are returned, even if the other is on the line of purple.
+    The first point is always in the dominant direction. If the dominant cannot
+    be found, the complementary will be used and indicated with a negative sign.
+
+    If `reverse` is set, then the complementary wavelength is returned instead,
+    with the points arranged favoring the complementary wavelength. If it cannot
+    be found, the dominant wavelength is returned with a negative sign.
+
+    If `closet` is set, wavelengths are rounded to the closest.
+    """
+
+    w1 = w2 = math.nan
+    dominant = [math.nan, math.nan]
+    complementary = [math.nan, math.nan]
+
+    # Achromatic, no wavelength
+    if all(abs(a - b) < 1e-12 for a, b in zip(xy, white)):
+        return w1, dominant, complementary
+
+    # Look for first intersection of the line drawn through the white point
+    # and the current color with the spectral locus. Check the dominant and
+    # complementary, but return as soon as we have the dominant. If no dominant
+    # is found, we'll use the complementary.
+    locus = [util.xyz_to_xyY(cmfs.CIE_1931_2DEG[r], white)[:-1] for r in range(LOCUS_START, LOCUS_END + 1, LOCUS_STEP)]
+    start = xy_to_angle(locus[0], white)
+    current = xy_to_angle(xy, white, start)
+    invert = xy_to_angle(xy, white, start, invert=True)
+    found = [False, False]
+    for i in range(1, len(locus) - 2):
+        # Get the next locus point angle
+        a_next = xy_to_angle(locus[i], white, start)
+
+        # Check if our angle is greater than the current locus point's angle
+        for j in range(0, 2):
+
+            # If has already been found, skip
+            target = invert if j else current
+            if a_next > target or found[j]:
+                continue
+
+            # Previous index
+            i0 = i - 1
+
+            # Get the intersection
+            if target == a_next:
+                intersect = locus[i]  # type: Vector | None
+            elif target == xy_to_angle(locus[i0], white, start):
+                intersect = locus[i0]
+            else:
+                intersect = alg.ray_line_intersect(white, xy, locus[i0], locus[i])
+            if intersect is None:  # pragma: no cover
+                continue
+
+            # Use the intersection to estimate the interpolation factor for the wavelength,
+            # and then get the interpolated value via Sprague interpolation.
+            i2 = 0 if abs(locus[i0][0] - locus[i][0]) > abs(locus[i0][1] - locus[i][1]) else 1
+            f = alg.ilerp(locus[i0][i2], locus[i][i2], intersect[i2])
+            w = alg.lerp(LOCUS_START + i0, LOCUS_START + i, f)
+            intersect = util.xyz_to_xyY(cmfs.CIE_1931_2DEG[w], white)[:-1]
+
+            if j == 0:
+                dominant = intersect
+                w1 = w
+                found[j] = True
+                if not reverse:
+                    break
+            else:
+                complementary = intersect
+                w2 = w
+                found[j] = True
+                if reverse:
+                    break
+
+        if found[reverse]:
+            break
+
+    # Unlikely catastrophic failure
+    if not any(found):  # pragma: no cover
+        return w1, dominant, complementary
+
+    # Swap dominant and complementary if we are looking for complementary
+    if reverse:
+        dominant, complementary = complementary, dominant
+        w1, w2 = w2, w1
+
+    # If dominant isn't found, it is on the line of purples; use complementary instead
+    if not found[reverse]:
+        pt = alg.ray_line_intersect(white, xy, locus[0], locus[-1])
+        # Shouldn't happen, but just in case
+        if pt is not None:  # pragma: no cover
+            dominant = pt
+        w1 = -alg.round_half_up(w2) if closest else -w2
+    else:
+        complementary = dominant
+        if closest:
+            w1 = alg.round_half_up(w1)
+
+    return w1, dominant, complementary
+
+
+def wavelength_to_color(wavelength: float) -> Vector:
+    """Return the XYZ value for the specified wavelength."""
+
+    if wavelength < LOCUS_START or wavelength > LOCUS_END:
+        raise ValueError(f'{wavelength}nm exceeds the range of {LOCUS_START}nm - {LOCUS_END}nm')
+
+    # Wavelength is within the CMFs
+    return cmfs.CIE_1931_2DEG[wavelength]

docs/src/dictionary/en-custom.txt

@@ -395,6 +395,7 @@ undersaturated
 uniformities
 unintuitive
 unsaturate
+unscaled
 uv
 uvY
 validator

docs/src/markdown/.snippets/abbr.md

@@ -19,4 +19,4 @@
 *[SDR]: Standard Dynamic Range
 *[SVG]: Scalable Vector Graphics
 *[UCS]: Uniform Chromaticity Space
-*[WCAG]: Web Content Accessibility Guidelines
+*[WCAG]: Web Content Accessibility Guidelines

docs/src/markdown/images/dominant-wavelength-cyan.png

@@ -90,8 +90,7 @@ out to even this limit.
 ## Out of Gamut Temperatures
 
 It should be noted that `blackbody()` normalizes/scales the returned colors by default as the colors are often much too
-bright initially, all having a max luminance. This scaling is usually done under a linear RGB color space, Linear
-Rec2020 being the default as it encompasses the entire curve.
+bright initially, all having a max luminance. This scaling is done is linear sRGB.
 
 Keep in mind that if the color is not in the display gamut it will need to be gamut mapped, and the gamut mapped value
 will not exhibit the same temperature. How far off it is will be depends on the disparity of the gamut sizes and how the