Add a luminance preserving version of scale GMA (#506)

Author: facelessuser

coloraide/__meta__.py

@@ -204,5 +204,5 @@ def parse_version(ver: str) -> Version:
     return Version(major, minor, micro, release, pre, post, dev)
 
 
-__version_info__ = Version(8, 4, 0, "final")
+__version_info__ = Version(8, 5, 0, "final")
 __version__ = __version_info__._get_canonical()

coloraide/algebra.py

@@ -1182,58 +1182,90 @@ def pprint(value: float | ArrayLike) -> None:
     print(pretty(value))
 
 
-def ray_line_intersect(
-    a1: VectorLike,
-    a2: VectorLike,
-    b1: VectorLike,
-    b2: VectorLike,
+def point_on_segment(a: VectorLike, b: VectorLike, p: VectorLike, abs_tol: float = ATOL) -> bool:
+    """Point on line segment."""
+
+    l = len(p)
+
+    ab = [b[i] - a[i] for i in range(l)]
+    ap = [p[i] - a[i] for i in range(l)]
+
+    # Check if `ap` is parallel to `ab` (cross product is zero)
+    cp = cross(ab, ap)
+    if _any(abs(c) > abs_tol for c in cp):
+        return False
+
+    # Check if the point is between a and b
+    for i in range(l):
+        if abs(ab[i]) < abs_tol:
+            continue
+        t = ap[i] / ab[i]
+        break
+
+    # See if `ab` is a point and `p` is that point
+    else:
+        return _all(abs(i - j) < abs_tol for i, j in zip(a, p))
+
+    return 0 <= t <= 1
+
+
+def line_interesect(
+    s1: VectorLike,
+    e1: VectorLike,
+    s2: VectorLike,
+    e2: VectorLike,
+    rel_tol: float = RTOL,
     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
+    """
+    Find intersection of two lines.
+
+    This was designed particularly for 3D intersection, but can be used for either 2D or 3D,
+    but 2D line intersection could be calculated with less work using other methods if performance
+    was of importance.
+
+    3D lines rarely intersect, but often the shortest line between can be found.
+    If the shortest line is has no length (a point) then it is an actual intersection.
+    Our cases are constructed such that an intersection is expected, and a line is not sufficient.
+    We can verify closeness of the points (to account for floating point errors) to verify that within
+    some expected threshold, the two line points are essentially a point and an intersection is found.
+    """
+
+    # Line segment difference
+    l1 = [a - b for a, b in zip(e1, s1)]
+    l2 = [a - b for a, b in zip(e2, s2)]
+
+    # Magnitude
+    m1 = math.sqrt(sum(a ** 2 for a in l1))
+    m2 = math.sqrt(sum(a ** 2 for a in l2))
+
+    # One of the lines is a point.
+    if m1 < abs_tol:
+        return list(s1) if point_on_segment(e2, s2, s1, abs_tol=abs_tol) else None
+    elif m2 < abs_tol:
+        return list(s2) if point_on_segment(e1, s1, s2, abs_tol=abs_tol) else None
+
+    # Unit vector
+    u1 = [a / m1 for a in l1]
+    u2 = [a / m2 for a in l2]
+    # Direction projection
+    u = vdot(u1, u2)
+    if u == 1:  # pragma: no cover
         return None
-    t = dot(dap, dp, dims=D1) / denom
-    # Check if intersection is within bounds
-    if 0 <= t <= 1:
-        # Intersect
-        return add(b1, multiply(t, db, dims=SC_D1), dims=D1)
-    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
+    # Separation projection
+    sp = [a - b for a, b in zip(s2, s1)]
+    sp1 = vdot(sp, u1)
+    sp2 = vdot(sp, u2)
+    # Distance along lines
+    d1 = (sp1 - u * sp2) / (1 - u * u)
+    d2 = (sp2 - u * sp1) / (u * u - 1)
+    # Calculate points of closest line
+    p1 = [a + b for a, b in zip(s1, [x * d1 for x in u1])]
+    p2 = [a + b for a, b in zip(s2, [x * d2 for x in u2])]
+    # If points are close enough, assume intersect, otherwise raise error
+    if not _all(math.isclose(i, j, rel_tol=rel_tol, abs_tol=abs_tol) for i, j in zip(p1, p2)):  # 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
+    return p1
 
 
 def all(a: float | ArrayLike) -> bool:  # noqa: A001

coloraide/color.py

@@ -70,6 +70,7 @@
 from .gamut.fit_oklch_chroma import OkLChChroma
 from .gamut.fit_raytrace import RayTrace
 from .gamut.fit_scale import Scale
+from .gamut.fit_scale_luminance import ScaleLuminance
 from .cat import CAT, Bradford
 from .filters import Filter
 from .filters.w3c_filter_effects import Sepia, Brightness, Contrast, Saturate, Opacity, HueRotate, Grayscale, Invert
@@ -529,6 +530,7 @@ def blackbody(
         scale_space: str | None = None,
         max_saturation: bool = True,
         clip_negative: bool = False,
+        preserve_luminance: bool = False,
         **kwargs: Any
     ) -> Self:
         """
@@ -553,7 +555,8 @@ def blackbody(
             scale=scale,
             scale_space=scale_space,
             max_saturation=max_saturation,
-            clip_negative=clip_negative
+            clip_negative=clip_negative,
+            preserve_luminance=preserve_luminance
         )
 
     def cct(self, *, method: str | None = None, **kwargs: Any) -> Vector:
@@ -844,7 +847,8 @@ def chromaticity(
         scale: bool = False,
         scale_space: str | None = None,
         max_saturation: bool = False,
-        clip_negative: bool = False
+        clip_negative: bool = False,
+        preserve_luminance: bool = False
     ) -> Self:
         """
         Create a color from chromaticity coordinates.
@@ -882,7 +886,8 @@ def chromaticity(
                 color,
                 scale_space=scale_space if scale_space is not None else 'srgb-linear',
                 max_saturation=max_saturation,
-                clip_negative=clip_negative
+                clip_negative=clip_negative,
+                preserve_luminance=preserve_luminance
             )
 
         # Convert to targeted color space
@@ -1382,7 +1387,8 @@ def from_wavelength(
         scale: bool = True,
         scale_space: str | None = None,
         max_saturation: bool = True,
-        clip_negative: bool = False
+        clip_negative: bool = False,
+        preserve_luminance: bool = False
     ) -> Self:
         """Create a color from a wavelength."""
 
@@ -1394,7 +1400,8 @@ def from_wavelength(
             scale=scale,
             scale_space=scale_space,
             max_saturation=max_saturation,
-            clip_negative=clip_negative
+            clip_negative=clip_negative,
+            preserve_luminance=preserve_luminance
         )
 
     @overload
@@ -1645,6 +1652,7 @@ def alpha(
         OkLChChroma(),
         RayTrace(),
         Scale(),
+        ScaleLuminance(),
 
         # Filters
         Sepia(),

coloraide/gamut/__init__.py

@@ -112,25 +112,46 @@ def from_base(self, coords: Vector) -> Vector:
     return RGB()
 
 
-def adjust_luminance(color: Color, Y: float, white: VectorLike) -> None:
+def adjust_luminance(
+    color: Color,
+    Y: float,
+    white: VectorLike,
+    max_luminance: float = 1.0,
+    preserve_luminance: bool = True
+) -> None:
     """Adjust luminance of a color."""
 
     with color.within('xyz-d65') as c:
-        c.convert('xyz-d65', in_place=True)
         d65 = c._space.WHITE
         adapt = d65 != white
         xyz = c.chromatic_adaptation(d65, white, c[:-1]) if adapt else c[:-1]
+        Y = alg.clamp(Y, 0.0, max_luminance)
+        # Luminance below the cusp can just be restored
         if xyz[1] > Y:
             xyz = util.xy_to_xyz(util.xyz_to_xyY(xyz, white)[:-1], Y)
             c[:-1] = c.chromatic_adaptation(white, d65, xyz) if adapt else xyz
+        # Luminance above the cusp requires us to find the intersection of the vectors of the
+        # path between the color and white and those same colors with the adjusted luminance.
+        elif preserve_luminance and xyz[1] < Y:
+            xyy = util.xyz_to_xyY(xyz, white)
+            intersect = alg.line_interesect(
+                xyz,
+                util.xy_to_xyz(white, max_luminance),
+                util.xy_to_xyz(xyy[:2], Y),
+                util.xy_to_xyz(white, Y)
+            )
+            # Update color if we found an intersection
+            if intersect is not None:
+                c[:-1] = c.chromatic_adaptation(white, d65, intersect) if adapt else intersect
 
 
 def scale_rgb(
     color: Color,
     *,
     scale_space: str,
     clip_negative: bool = False,
-    max_saturation: bool = False
+    max_saturation: bool = False,
+    preserve_luminance: bool = False
 ) -> None:
     """Apply color scaling."""
 
@@ -171,10 +192,9 @@ def scale_rgb(
         rgb[i] = alg.clamp((rgb[i] - mn) / mx if mx else (rgb[i] - mn), 0.0, 1.0) * maximum
     mapcolor[:-1] = cs.to_base(rgb) if coerced else rgb
 
-    # If the current luminance is greater than the original luminance,
-    # set the luminance to the original. Set it xyY with the same white point.
+    # Check if luminance doesn't match and update accordingly.
     if not max_saturation:
-        adjust_luminance(mapcolor, Y, white)
+        adjust_luminance(mapcolor, Y, white, maximum, preserve_luminance)
 
     # Clip in the original gamut bound color space and update the original color
     clip_channels(mapcolor.convert(orig_space, in_place=True))

coloraide/gamut/fit_scale_luminance.py

@@ -0,0 +1,25 @@
+"""Luminance preserving scale."""
+from __future__ import annotations
+from .fit_scale import Scale
+from . import scale_rgb
+from typing import TYPE_CHECKING, Any
+
+if TYPE_CHECKING:  # pragma: no cover
+    from .. import Color
+
+
+class ScaleLuminance(Scale):
+    """Luminance preserving scale."""
+
+    NAME = 'scale-luminance'
+
+    def fit(
+        self,
+        color: Color,
+        space: str,
+        preserve_luminance: bool = True,
+        **kwargs: Any
+    ) -> None:
+        """Scale the color within its gamut but preserve L and h as much as possible."""
+
+        scale_rgb(color, scale_space=space, preserve_luminance=preserve_luminance, **kwargs)

coloraide/spectrum.py

@@ -38,6 +38,31 @@ def xy_to_angle(xy: VectorLike, white: VectorLike, offset: float = 0.0, invert:
     return angle
 
 
+def ray_line_intersect(
+    a1: VectorLike,
+    a2: VectorLike,
+    b1: VectorLike,
+    b2: VectorLike,
+    abs_tol: float = alg.ATOL
+) -> Vector | None:
+    """Find the intersection of a 2D ray and 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 = alg.dot(dap, db, dims=alg.D1)
+    # Parallel cases
+    if abs(denom) < abs_tol:  # pragma: no cover
+        return None
+    t = alg.dot(dap, dp, dims=alg.D1) / denom
+    # Check if intersection is within bounds
+    if 0 <= t <= 1:
+        # Intersect
+        return alg.add(b1, alg.multiply(t, db, dims=alg.SC_D1), dims=alg.D1)
+    return None  # pragma: no cover
+
+
 def closest_wavelength(
     xy: VectorLike,
     white: VectorLike = WHITE,
@@ -103,7 +128,7 @@ def closest_wavelength(
             elif target == xy_to_angle(locus[i0], white, start):
                 intersect = locus[i0]
             else:
-                intersect = alg.ray_line_intersect(white, xy, locus[i0], locus[i])
+                intersect = ray_line_intersect(white, xy, locus[i0], locus[i])
             if intersect is None:  # pragma: no cover
                 continue
 
@@ -141,7 +166,7 @@ def closest_wavelength(
 
     # 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])
+        pt = 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

docs/src/markdown/about/changelog.md

@@ -3,6 +3,11 @@ icon: lucide/scroll-text
 ---
 # Changelog
 
+## 8.5
+
+-   **NEW**: Add `perserve_luminance` option to `scale` to allow for a luminance preserving version.
+-   **NEW**: Add `scale-luminance` gamut mapping method as a pre-configured luminance preserving scale method.
+
 ## 8.4
 
 -   **NEW**: Add a new gamut mapping approach called `scale` that scales colors within a linear RGB space and can