diff --git a/coloraide/gamut/fit_raytrace.py b/coloraide/gamut/fit_raytrace.py
index 0e68f3c3..f6ffec72 100644
--- a/coloraide/gamut/fit_raytrace.py
+++ b/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:
diff --git a/docs/src/markdown/about/changelog.md b/docs/src/markdown/about/changelog.md
index c9ed3d06..86b93dd1 100644
--- a/docs/src/markdown/about/changelog.md
+++ b/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
diff --git a/docs/src/markdown/gamut.md b/docs/src/markdown/gamut.md
index 6f1fc21e..75de15c3 100644
--- a/docs/src/markdown/gamut.md
+++ b/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
diff --git a/docs/src/markdown/images/raytrace-gma-improve.png b/docs/src/markdown/images/raytrace-gma-improve.png
new file mode 100644
index 00000000..e840252e
Binary files /dev/null and b/docs/src/markdown/images/raytrace-gma-improve.png differ
diff --git a/tests/test_a98_rgb.py b/tests/test_a98_rgb.py
index 2c55dfab..6e53a4e5 100644
--- a/tests/test_a98_rgb.py
+++ b/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)')
     ]
diff --git a/tests/test_display_p3.py b/tests/test_display_p3.py
index c673bdfb..b979a87b 100644
--- a/tests/test_display_p3.py
+++ b/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)')
     ]
diff --git a/tests/test_gamut.py b/tests/test_gamut.py
index 2a4d4be2..408f90c7 100644
--- a/tests/test_gamut.py
+++ b/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."""
diff --git a/tests/test_prophoto_rgb.py b/tests/test_prophoto_rgb.py
index e946fa4e..7d3041e7 100644
--- a/tests/test_prophoto_rgb.py
+++ b/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)')
     ]
diff --git a/tests/test_rec2020.py b/tests/test_rec2020.py
index f1c9ce5b..be9df5bd 100644
--- a/tests/test_rec2020.py
+++ b/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)')
     ]
diff --git a/tests/test_rec709.py b/tests/test_rec709.py
index af0e24dd..8a5e779d 100644
--- a/tests/test_rec709.py
+++ b/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)')
     ]
diff --git a/tests/test_srgb.py b/tests/test_srgb.py
index f3c35499..6de7c8df 100644
--- a/tests/test_srgb.py
+++ b/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)')
     ]
diff --git a/tools/gamutcheck.py b/tools/gamutcheck.py
index 24e3a683..90e2cc58 100644
--- a/tools/gamutcheck.py
+++ b/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
 
diff --git a/tools/raytrace_test.py b/tools/raytrace_test.py
index 339367be..1d65e664 100644
--- a/tools/raytrace_test.py
+++ b/tools/raytrace_test.py
@@ -2,6 +2,7 @@
 import sys
 import os
 import plotly.graph_objects as go
+import plotly.io as io
 import argparse
 import json
 
@@ -262,6 +263,10 @@ def simulate_raytrace_gamut_mapping(args):
     else:
         gamutcolor = mapcolor.convert(space)
 
+        # Threshold for anchor adjustment
+        low = 1e-6
+        high = bmx - low
+
         # Create a ray from our current color to the color with zero chroma.
         # Trace the line to the RGB cube finding the intersection.
         # In between iterations, correct the L and H and then cast a ray
@@ -297,7 +302,13 @@ def simulate_raytrace_gamut_mapping(args):
                         )
 
                 gamutcolor.convert(space, in_place=True)
-            intersection = raytrace_box(achromatic, gamutcolor[:-1], bmax=bmax)
+
+            coords = gamutcolor[:-1]
+            intersection = raytrace_box(achromatic, coords, bmax=bmax)
+
+            if i and all(low < x < high for x in coords):
+                achromatic = gamutcolor[:-1]
+
             if intersection:
                 points.append(gamutcolor[:-1])
                 points.append(intersection)
@@ -403,6 +414,7 @@ def simulate_raytrace_gamut_mapping(args):
     )
     parser.add_argument('--height', '-H', type=int, default=800, help="Diagram height.")
     parser.add_argument('--width', '-W', type=int, default=800, help="Diagram width.")
+    parser.add_argument('--output', '-o', default='', help='Output file.')
     args = parser.parse_args()
 
     if args.gamut_color:
@@ -410,4 +422,18 @@ def simulate_raytrace_gamut_mapping(args):
     else:
         fig = raytrace(args)
 
-    fig.show()
+
+    # Show or save the data as an image, etc.
+    if fig:
+        if args.output:
+            filetype = os.path.splitext(args.output)[1].lstrip('.').lower()
+            if filetype == 'html':
+                with open(args.output, 'w') as f:
+                    f.write(io.to_html(fig))
+            elif filetype == 'json':
+                io.write_json(fig, args.output)
+            else:
+                with open(args.output, 'wb') as f:
+                    f.write(fig.to_image(format=filetype))
+        else:
+            fig.show()