Improve the ray trace algorithm to yield better results in all spaces (…

…#435)

The ray trace approach was designed focusing mainly on Oklab/OkLCh.
Further evaluation shows that some spaces twist a little different
through the RGB spaces and can yield results not as good, CIELCh
being an example.

Adjust the algorithm so that if a point is found along the reduction
path below the gamut surface, that point becomes the new anchor.

This generally tightness the range in a few passes and prevent some
spaces from having larger than expected hue deviations.

coloraide/gamut/fit_raytrace.py

@@ -20,25 +20,22 @@
     from ..color import Color
 
 
-def project_onto(p1: Vector, p2: Vector, p0: Vector) -> Vector:
+def project_onto(a: Vector, b: Vector, o: Vector) -> Vector:
     """
     Using 3 points, create two vectors with a shared origin and project the first vector onto the second.
 
-    - `p1`:  point used to define the magnitude of the first vector (`v1`) with origin `p0`.
-    - `p2`:  point used to define the magnitude of the second vector (`v2`) with origin `p0`.
-    - `p0`:  the origin point of both `v1` and `v2`.
+    - `a`:  point used to define the head of the first vector `OA`.
+    - `b`:  point used to define the head of the second vector `OB`.
+    - `o`:  the origin/tail point of both vector `OA` and `OB`.
     """
 
-    # Separate into points
-    x1, y1, z1 = p1
-    x2, y2, z2 = p2
-    x3, y3, z3 = p0
     # Create vector from points
-    v1 = [x1 - x3, y1 - y3, z1 - z3]
-    v2 = [x2 - x3, y2 - y3, z2 - z3]
-    # Project v1 onto v2 and convert back to a point
-    r = alg.vdot(v1, v2) / alg.vdot(v2, v2)
-    return [v2[0] * r + x3, v2[1] * r + y3, v2[2] * r + z3]
+    ox, oy, oz = o
+    vec_oa = [a[0] - ox, a[1] - oy, a[2] - oz]
+    vec_ob = [b[0] - ox, b[1] - oy, b[2] - oz]
+    # Project `vec_oa` onto `vec_ob` and convert back to a point
+    r = alg.vdot(vec_oa, vec_ob) / alg.vdot(vec_ob, vec_ob)
+    return [vec_ob[0] * r + ox, vec_ob[1] * r + oy, vec_ob[2] * r + oz]
 
 
 def hwb_to_srgb(coords: Vector) -> Vector:  # pragma: no cover
@@ -282,25 +279,37 @@ def fit(
             # In between iterations, correct the L and H and then cast a ray
             # to the new corrected color finding the intersection again.
             mapcolor.convert(space, in_place=True)
+
+            # Interpolation path
+            if adaptive:
+                start = [light, *ab]
+                end = [alight, 0.0, 0.0]
+
+            # Threshold for anchor adjustment
+            low = 1e-6
+            high = bmx - low
+
             for i in range(4):
                 if i:
                     mapcolor.convert(pspace, in_place=True, norm=False)
 
+                    # Correct the point onto the desired interpolation path
                     if adaptive:
-                        # Correct the point onto the desired interpolation path
                         if polar:
                             mapcolor[l], a_, b_ = project_onto(
                                 [mapcolor[l], *alg.polar_to_rect(mapcolor[c], mapcolor[h])],
-                                [light, *ab],
-                                [alight, 0.0, 0.0]
+                                start,
+                                end
                             )
                             mapcolor[c], mapcolor[h] = alg.rect_to_polar(a_,b_)
                         else:
                             mapcolor[l], mapcolor[a], mapcolor[b] = project_onto(
                                 [mapcolor[l], mapcolor[a], mapcolor[b]],
-                                [light, *ab],
-                                [alight, 0.0, 0.0]
+                                start,
+                                end
                             )
+
+                    # Simple correction for constant lightness
                     else:
                         # Correct lightness and hue
                         mapcolor[l] = alight
@@ -314,8 +323,14 @@ def fit(
 
                     mapcolor.convert(space, in_place=True)
 
-                # Cast a ray to our anchor point.
-                intersection = raytrace_box(anchor, mapcolor[:-1], bmax=bmax)
+                coords = mapcolor[:-1]
+                intersection = raytrace_box(anchor, coords, bmax=bmax)
+
+                # Adjust anchor point closer to surface to improve results for some spaces.
+                # Don't move point too close to the surface to avoid corner cases with some spaces.
+                # OkLCh/Oklab does not require this.
+                if i and all(low < x < high for x in coords):
+                    anchor = mapcolor[:-1]
 
                 # Update color with the intersection point on the RGB surface.
                 if intersection:

docs/src/markdown/about/changelog.md

@@ -29,6 +29,7 @@
     methods such as `interpolate()`, etc.
 -   **NEW**: `compose()` has been deprecated in favor of the new `layer()` method and will be removed at some future
     time.
+-   **NEW**: Improve experimental `raytrace` gamut mapping approach when performed in certain perceptual spaces.
 -   **BREAK**: The experimental `raytrace` gamut mapping method now uses OkLCh by default instead of CIELCh (D65).
 -   **BREAK**: Pre-configured `oklch-raytrace` and `lch-raytrace` variants of the experimental `raytrace` gamut mapping
     method have been removed to reduce included plugins. OkLCh is the default now and users can still specify CIELCh and

docs/src/markdown/gamut.md

@@ -461,10 +461,11 @@ version will be used automatically for best results. Currently ColorAide can gam
 spaces as they either have an RGB gamut or can be coerced into one.
 
 The ray trace approach works by taking a given color and converting it to a perceptual Lab-ish or LCh-ish color space
-(the default being OkLCh) and then calculates the achromatic version of the color (chroma set to zero). Assuming the
-lightness is within bounds, a ray is cast from the inside of the cube, from the achromatic point to the current color.
-The intersection along this path with the RGB gamut surface is then found.  If the achromatic color exceeds the maximum
-or minimum lightness of the gamut, the respective maximum or minimum achromatic color is returned. 
+(the default being OkLCh) and then calculates the achromatic version of the color (chroma set to zero) which will be our
+anchor point. Assuming our anchor point is within bounds, a ray is cast from the inside of the cube, from the anchor
+point to the current color. The intersection along this path with the RGB gamut surface is then found.  If the
+achromatic color exceeds the maximum or minimum lightness of the gamut, the respective maximum or minimum achromatic
+color is returned. 
 
 /// note | Ray Trace Algorithm
 The ray trace algorithm is based on the [slab method](https://en.wikipedia.org/wiki/Slab_method). The intersection that
@@ -474,20 +475,26 @@ end point.
 
 The intersection of the line and the gamut surface represents an approximation of the the most saturated color for that
 lightness and hue, but because the RGB space is not perceptual, the initial approximation is likely to be off because
-decreasing chroma and holding lightness and hue constant in a perceptual space will create curved path through the
+decreasing chroma and holding lightness and hue constant in a perceptual space will create a curved path through the
 RGB space. In order to converge on a point as close as possible to the actual most saturated color with the given hue
 and lightness, we must refine our result with a few additional iterations.
 
 In order to converge on the actual chroma reduced color we seek, we can take the first intersection we find and correct
-the color in the perceptual color space by setting the hue and lightness back to the original color's hue and lightness.
-The corrected color becomes our new current color and should be a much closer color on the reduced chroma line. We can
-repeat this process a few more times, each time finding a better, closer color on the path. After about three
-_additional_ iterations (a combined total of four for the entire process), we will be close enough where we can stop.
-Finally, we can then clip off floating point math errors. With this, we will now have a more accurate approximation of
-the color we seek.
+the color in the perceptual color space by projecting the point back onto the chroma reduction path, correcting the
+color's hue and lightness. The corrected color becomes our new current color and should be a much closer color on the
+reduced chroma line. We can repeat this process a few more times, each time finding a better, closer color on the path.
+After about three _additional_ iterations (a combined total of four for the entire process), we will be close enough
+where we can stop. Finally, we can then clip off floating point math errors. With this, we will now have a more accurate
+approximation of the color we seek.
 
 ![Ray Trace Gamut Mapping Example](images/raytrace-gma.png)
 
+One final improvement is that during the correction step, where we adjust surface point back onto the chroma reduction
+path, if we find a point below the gamut surface, we can adjust our anchor to be this new point, closer to the gamut
+surface, which in some spaces will help to converge closer to our ideal color than they would without the adjustment.
+
+![Ray Trace Gamut Mapping Example](images/raytrace-gma-improve.png)
+
 /// note
 For accuracy, iterations could be increased further which would reduce a potential ∆h shift even more, but ColorAide has
 opted to keep iterations at 4 which can gamut map colors to sRGB with ∆h shift of less than 1, and when gamut mapping

tests/test_a98_rgb.py

@@ -54,7 +54,7 @@ class TestA98RGBSerialize(util.ColorAssertsPyTest):
         (
             'color(a98-rgb 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(a98-rgb 1 0.53016 0.38997)'
+            'color(a98-rgb 1 0.5301 0.3906)'
         ),
         ('color(a98-rgb 1.2 0.2 0)', {'fit': False}, 'color(a98-rgb 1.2 0.2 0)')
     ]

tests/test_display_p3.py

@@ -54,7 +54,7 @@ class TestDisplayP3Serialize(util.ColorAssertsPyTest):
         (
             'color(display-p3 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(display-p3 1 0.48317 0.24021)'
+            'color(display-p3 1 0.48195 0.26354)'
         ),
         ('color(display-p3 1.2 0.2 0)', {'fit': False}, 'color(display-p3 1.2 0.2 0)')
     ]

tests/test_gamut.py

@@ -278,7 +278,7 @@ def test_raytrace_adaptive_lightness_lab(self):
 
         self.assertColorEqual(
             Color('display-p3', [1, 1, 0]).fit('srgb', **options),
-            Color('color(display-p3 0.89593 0.90035 0.29413)')
+            Color('color(display-p3 0.89593 0.90035 0.29412)')
         )
 
     def test_raytrace_adaptive_lightness_lch(self):
@@ -295,7 +295,7 @@ def test_raytrace_adaptive_lightness_lch(self):
 
         self.assertColorEqual(
             Color('display-p3', [1, 1, 0]).fit('srgb', **options),
-            Color('color(display-p3 0.89593 0.90035 0.29413)')
+            Color('color(display-p3 0.89593 0.90035 0.29412)')
         )
 
     def test_sdr_extremes_low(self):
@@ -341,7 +341,7 @@ def test_trace_face2(self):
         """Test tracing of face 2."""
 
         color = Color('oklch(90% .4 60)')
-        self.assertColorEqual(color.fit('srgb', method='raytrace', pspace='oklch'), Color('oklch(0.9 0.06478 60.008)'))
+        self.assertColorEqual(color.fit('srgb', method='raytrace', pspace='oklch'), Color('oklch(0.9 0.06478 60)'))
 
     def test_trace_face3(self):
         """Test tracing of face 3."""
@@ -359,7 +359,7 @@ def test_trace_face5(self):
         """Test tracing of face 5."""
 
         color = Color('oklch(30% .4 100)')
-        self.assertColorEqual(color.fit('srgb', method='raytrace', pspace='oklch'), Color('oklch(0.3 0.06237 99.999)'))
+        self.assertColorEqual(color.fit('srgb', method='raytrace', pspace='oklch'), Color('oklch(0.3 0.06237 100)'))
 
     def test_trace_face6(self):
         """Test tracing of face 6."""

tests/test_prophoto_rgb.py

@@ -54,7 +54,7 @@ class TestProPhotoRGBSerialize(util.ColorAssertsPyTest):
         (
             'color(prophoto-rgb 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(prophoto-rgb 1 0.40937 0.1277)'
+            'color(prophoto-rgb 1 0.40879 0.18246)'
         ),
         ('color(prophoto-rgb 1.2 0.2 0)', {'fit': False}, 'color(prophoto-rgb 1.2 0.2 0)')
     ]

tests/test_rec2020.py

@@ -54,7 +54,7 @@ class TestRec2020Serialize(util.ColorAssertsPyTest):
         (
             'color(rec2020 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(rec2020 1 0.46393 0.16096)'
+            'color(rec2020 1 0.46248 0.19988)'
         ),
         ('color(rec2020 1.2 0.2 0)', {'fit': False}, 'color(rec2020 1.2 0.2 0)')
     ]

tests/test_rec709.py

@@ -54,7 +54,7 @@ class TestRec709Serialize(util.ColorAssertsPyTest):
         (
             'color(--rec709 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(--rec709 1 0.41123 0.24137)'
+            'color(--rec709 1 0.41112 0.24303)'
         ),
         ('color(--rec709 1.2 0.2 0)', {'fit': False}, 'color(--rec709 1.2 0.2 0)')
     ]

tests/test_srgb.py

@@ -129,7 +129,7 @@ class TestsRGBSerialize(util.ColorAssertsPyTest):
         (
             'color(srgb 1.2 0.2 0)',
             {'color': True, 'fit': {'method': 'raytrace', 'pspace': 'lch-d65'}},
-            'color(srgb 1 0.4597 0.31491)'
+            'color(srgb 1 0.45963 0.31602)'
         ),
         ('color(srgb 1.2 0.2 0)', {'color': True, 'fit': False}, 'color(srgb 1.2 0.2 0)')
     ]

tools/gamutcheck.py

@@ -9,11 +9,12 @@
 from coloraide.everything import ColorAll as Color  # noqa: E402
 
 
-def run(gamut, space, max_chroma, projection, adaptive, calc_near):
+def run(gamut, space, max_chroma, adaptive, calc_near):
     """Test the gamut mapping algorithm noting the ∆L, ∆h, and the worst ∆h offender."""
 
     far_delta_h = max_delta_h = 0
     far_delta_l = max_delta_l = 0
+    far_de = max_de = 0
     far_worst = worst = None
     lch = Color('white').convert(space)
     l, _, h = lch._space.indexes()
@@ -24,11 +25,15 @@ def run(gamut, space, max_chroma, projection, adaptive, calc_near):
         end = light
         for hue in range(360):
             c1 = Color(space, [light * scale, max_chroma, hue / 2])
-            c2 = c1.clone().fit(gamut, method='raytrace', pspace=space, projection=projection, adaptive=adaptive)
+            c2 = c1.clone().fit(gamut, method='raytrace', pspace=space, adaptive=adaptive)
+            de = c1.delta_e(c2, method='2000')
             dl = (c1[l] - c2[l])
             if abs(dl) > abs(max_delta_l):
                 max_delta_l = dl
 
+            if de > max_de:
+                max_de = de
+
             h1, h2 = c1[h], c2[h]
             dh = (h1 - h2)
             if dh > 180:
@@ -39,25 +44,31 @@ def run(gamut, space, max_chroma, projection, adaptive, calc_near):
                 max_delta_h = dh
                 worst = c1
 
-    print(f'=== Far: Chroma {max_chroma} (Lightness {start} - {end}) ===')
-    print('∆L =', max_delta_l)
-    print('∆h =', max_delta_h)
-    print('Worst ∆h offender =>', worst)
-
-    if abs(far_delta_h) > abs(max_delta_h):
-        max_delta_h = far_delta_h
-        worst = far_worst
-    if abs(far_delta_l) > abs(max_delta_l):
-        max_delta_l = far_delta_l
-
-    far_delta_h = max_delta_h
-    far_delta_l = max_delta_l
-    far_worst = worst
-
-    max_delta_h = 0
-    max_delta_l = 0
-    worst = None
-    start = end + 1
+        if light in (25, 50, 75, 99):
+            print(f'=== Far: Chroma {max_chroma} (Lightness {start} - {end}) ===')
+            print('∆L =', max_delta_l)
+            print('∆h =', max_delta_h)
+            print('∆E =', max_de)
+            print('Worst ∆h offender =>', worst)
+
+            if abs(far_delta_h) > abs(max_delta_h):
+                max_delta_h = far_delta_h
+                worst = far_worst
+            if abs(far_delta_l) > abs(max_delta_l):
+                max_delta_l = far_delta_l
+            if far_de > max_de:
+                max_de = far_de
+
+            far_delta_h = max_delta_h
+            far_delta_l = max_delta_l
+            far_de = max_de
+            far_worst = worst
+
+            max_delta_h = 0
+            max_delta_l = 0
+            max_de = 0
+            worst = None
+            start = end + 1
 
     if calc_near:
         for r, g, b in it.product(range(101), range(101), range(101)):
@@ -112,17 +123,14 @@ def main():
         '--max-chroma', '-c', type=float, default=0.8, help="Max chroma to test GMA with."
     )
     parser.add_argument(
-        '--projection', '-p', default='constant', help="Projection mode for the lightness."
-    )
-    parser.add_argument(
-        '--adaptive', '-a', type=float, default=0.05, help="Adaptive value."
+        '--adaptive', '-a', type=float, default=0.0, help="Adaptive value."
     )
     parser.add_argument(
         '--near', '-n', action='store_true', help="Calculate deltas when out of gamut color is close to boundary."
     )
     args = parser.parse_args()
 
-    run(args.gamut, args.lch, args.max_chroma, args.projection, args.adaptive, args.near)
+    run(args.gamut, args.lch, args.max_chroma, args.adaptive, args.near)
 
     return 0