How to find a maximum chroma value in the OKLCH color space for given hue and lightness?

Question

Chroma value varies for different hues, and I want to understand how I can find the maximum value that is displayable for the given hue and lightness programmatically.

I saw a couple of color pickers that can limit the selection visually, but I don't have access to their code, but anyway, I think it should be possible somehow.

Answer 1

Björn, who invented Oklab and Oklch, has published an article about gamut clipping using Oklch. In there you will find code that finds the bounds of sRGB in Oklch. The function gamut_clip_preserve_chroma should do pretty much what you want, so long as you strip off the RGB conversion steps at the beginning and end, and pretend that you have a ridiculously large C value — 1.0 will do.

I happen to need the same thing (well, a very similar thing: I want to make an RGB color as saturated as possible), so here's the crude TypeScript translation (yes it works):

// This is just https://bottosson.github.io/posts/gamutclipping/ crudely translated to TS

type Lab = { L: number; a: number; b: number };
type RGB = { r: number; g: number; b: number };
type sRGB = { r: number; g: number; b: number };

function linear_srgb_to_oklab(c: RGB): Lab {
  const l = 0.4122214708 * c.r + 0.5363325363 * c.g + 0.0514459929 * c.b;
  const m = 0.2119034982 * c.r + 0.6806995451 * c.g + 0.1073969566 * c.b;
  const s = 0.0883024619 * c.r + 0.2817188376 * c.g + 0.6299787005 * c.b;

  const l_ = Math.cbrt(l);
  const m_ = Math.cbrt(m);
  const s_ = Math.cbrt(s);

  return {
    L: 0.2104542553 * l_ + 0.7936177850 * m_ - 0.0040720468 * s_,
    a: 1.9779984951 * l_ - 2.4285922050 * m_ + 0.4505937099 * s_,
    b: 0.0259040371 * l_ + 0.7827717662 * m_ - 0.8086757660 * s_,
  };
}

function oklab_to_linear_srgb(c: Lab): RGB {
  const l_ = c.L + 0.3963377774 * c.a + 0.2158037573 * c.b;
  const m_ = c.L - 0.1055613458 * c.a - 0.0638541728 * c.b;
  const s_ = c.L - 0.0894841775 * c.a - 1.2914855480 * c.b;

  const l = l_ * l_ * l_;
  const m = m_ * m_ * m_;
  const s = s_ * s_ * s_;

  return {
    r: +4.0767416621 * l - 3.3077115913 * m + 0.2309699292 * s,
    g: -1.2684380046 * l + 2.6097574011 * m - 0.3413193965 * s,
    b: -0.0041960863 * l - 0.7034186147 * m + 1.7076147010 * s,
  };
}

function srgb_to_linear_rgb(c: sRGB): RGB {
    return {
        r: c.r <= 0.04045 ? c.r / 12.92 : Math.pow((c.r + 0.055) / 1.055, 2.4),
        g: c.g <= 0.04045 ? c.g / 12.92 : Math.pow((c.g + 0.055) / 1.055, 2.4),
        b: c.b <= 0.04045 ? c.b / 12.92 : Math.pow((c.b + 0.055) / 1.055, 2.4),
    };
}

function linear_rgb_to_srgb(c: RGB): sRGB {
    return {
        r: c.r <= 0.0031308 ? 12.92 * c.r : 1.055 * Math.pow(c.r, 1 / 2.4) - 0.055,
        g: c.g <= 0.0031308 ? 12.92 * c.g : 1.055 * Math.pow(c.g, 1 / 2.4) - 0.055,
        b: c.b <= 0.0031308 ? 12.92 * c.b : 1.055 * Math.pow(c.b, 1 / 2.4) - 0.055,
    };
}

function rgb_as_hex(c: sRGB): string {
    return "#" + [c.r, c.g, c.b].map(x => Math.round(x * 255).toString(16).padStart(2, "0")).join("");
}

function linear_rgb_as_hex(c: RGB): string {
    return rgb_as_hex(linear_rgb_to_srgb(c));
}

type LC = { L: number; C: number };
function compute_max_saturation(a: number, b: number): number {
  let k0, k1, k2, k3, k4, wl, wm, ws;

  if (-1.88170328 * a - 0.80936493 * b > 1) {
    k0 = +1.19086277;
    k1 = +1.76576728;
    k2 = +0.59662641;
    k3 = +0.75515197;
    k4 = +0.56771245;
    wl = +4.0767416621;
    wm = -3.3077115913;
    ws = +0.2309699292;
  } else if (1.81444104 * a - 1.19445276 * b > 1) {
    k0 = +0.73956515;
    k1 = -0.45954404;
    k2 = +0.08285427;
    k3 = +0.12541070;
    k4 = +0.14503204;
    wl = -1.2684380046;
    wm = +2.6097574011;
    ws = -0.3413193965;
  } else {
    k0 = +1.35733652;
    k1 = -0.00915799;
    k2 = -1.15130210;
    k3 = -0.50559606;
    k4 = +0.00692167;
    wl = -0.0041960863;
    wm = -0.7034186147;
    ws = +1.7076147010;
  }

  let S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;

  let k_l = +0.3963377774 * a + 0.2158037573 * b;
  let k_m = -0.1055613458 * a - 0.0638541728 * b;
  let k_s = -0.0894841775 * a - 1.2914855480 * b;

  {
    let l_ = 1 + S * k_l;
    let m_ = 1 + S * k_m;
    let s_ = 1 + S * k_s;

    let l = l_ * l_ * l_;
    let m = m_ * m_ * m_;
    let s = s_ * s_ * s_;

    let l_dS = 3 * k_l * l_ * l_;
    let m_dS = 3 * k_m * m_ * m_;
    let s_dS = 3 * k_s * s_ * s_;

    let l_dS2 = 6 * k_l * k_l * l_;
    let m_dS2 = 6 * k_m * k_m * m_;
    let s_dS2 = 6 * k_s * k_s * s_;

    let f = wl * l + wm * m + ws * s;
    let f1 = wl * l_dS + wm * m_dS + ws * s_dS;
    let f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2;

    S = S - f * f1 / (f1 * f1 - 0.5 * f * f2);
  }

  return S;
}

function find_cusp(a: number, b: number): LC {
  let S_cusp = compute_max_saturation(a, b);

  let rgb_at_max = oklab_to_linear_srgb({ L: 1, a: S_cusp * a, b: S_cusp * b });
  let L_cusp = Math.cbrt(1 / Math.max(rgb_at_max.r, rgb_at_max.g, rgb_at_max.b));
  let C_cusp = L_cusp * S_cusp;

  return { L: L_cusp, C: C_cusp };
}

function find_gamut_intersection(a: number, b: number, L1: number, C1: number, L0: number): number {
  let cusp = find_cusp(a, b);

  let t: number;
  if ((L1 - L0) * cusp.C - (cusp.L - L0) * C1 <= 0) {
    t = (cusp.C * L0) / (C1 * cusp.L + cusp.C * (L0 - L1));
  } else {
    t = (cusp.C * (L0 - 1)) / (C1 * (cusp.L - 1) + cusp.C * (L0 - L1));

    {
      let dL = L1 - L0;
      let dC = C1;

      let k_l = +0.3963377774 * a + 0.2158037573 * b;
      let k_m = -0.1055613458 * a - 0.0638541728 * b;
      let k_s = -0.0894841775 * a - 1.2914855480 * b;

      let l_dt = dL + dC * k_l;
      let m_dt = dL + dC * k_m;
      let s_dt = dL + dC * k_s;

      {
        let L = L0 * (1 - t) + t * L1;
        let C = t * C1;

        let l_ = L + C * k_l;
        let m_ = L + C * k_m;
        let s_ = L + C * k_s;

        let l = l_ * l_ * l_;
        let m = m_ * m_ * m_;
        let s = s_ * s_ * s_;

        let ldt = 3 * l_dt * l_ * l_;
        let mdt = 3 * m_dt * m_ * m_;
        let sdt = 3 * s_dt * s_ * s_;

        let ldt2 = 6 * l_dt * l_dt * l_;
        let mdt2 = 6 * m_dt * m_dt * m_;
        let sdt2 = 6 * s_dt * s_dt * s_;

        let r = 4.0767416621 * l - 3.3077115913 * m + 0.2309699292 * s - 1;
        let r1 = 4.0767416621 * ldt - 3.3077115913 * mdt + 0.2309699292 * sdt;
        let r2 = 4.0767416621 * ldt2 - 3.3077115913 * mdt2 + 0.2309699292 * sdt2;

        let u_r = r1 / (r1 * r1 - 0.5 * r * r2);
        let t_r = -r * u_r;

        let g = -1.2684380046 * l + 2.6097574011 * m - 0.3413193965 * s - 1;
        let g1 = -1.2684380046 * ldt + 2.6097574011 * mdt - 0.3413193965 * sdt;
        let g2 = -1.2684380046 * ldt2 + 2.6097574011 * mdt2 - 0.3413193965 * sdt2;

        let u_g = g1 / (g1 * g1 - 0.5 * g * g2);
        let t_g = -g * u_g;

        let b = -0.0041960863 * l - 0.7034186147 * m + 1.7076147010 * s - 1;
        let b1 = -0.0041960863 * ldt - 0.7034186147 * mdt + 1.7076147010 * sdt;
        let b2 = -0.0041960863 * ldt2 - 0.7034186147 * mdt2 + 1.7076147010 * sdt2;

        let u_b = b1 / (b1 * b1 - 0.5 * b * b2);
        let t_b = -b * u_b;

        t_r = u_r >= 0 ? t_r : Number.MAX_VALUE;
        t_g = u_g >= 0 ? t_g : Number.MAX_VALUE;
        t_b = u_b >= 0 ? t_b : Number.MAX_VALUE;

        t += Math.min(t_r, Math.min(t_g, t_b));
      }
    }
  }

  return t;
}

function clamp(x: number, min: number, max: number): number {
  return x < min ? min : x > max ? max : x;
}

function sgn(x: number): number {
  return +(0 < x) - +(x < 0);
}

function gamut_clip_preserve_chroma(rgb: RGB): RGB {
  if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0) {
    return rgb;
  }

  let lab = linear_srgb_to_oklab(rgb);

  let L = lab.L;
  let eps = 0.00001;
  let C = Math.max(eps, Math.sqrt(lab.a * lab.a + lab.b * lab.b));
  let a_ = lab.a / C;
  let b_ = lab.b / C;

  let L0 = clamp(L, 0, 1);

  let t = find_gamut_intersection(a_, b_, L, C, L0);
  let L_clipped = L0 * (1 - t) + t * L;
  let C_clipped = t * C;

  return oklab_to_linear_srgb({ L: L_clipped, a: C_clipped * a_, b: C_clipped * b_ });
}

function maximum_chroma_for_lh(L: number, h: number): number {
  let a = Math.cos(h);
  let b = Math.sin(h);

  let L0 = clamp(L, 0, 1);

  let t = find_gamut_intersection(a, b, L, 1, L0);

  return t;
}

function maximize_chroma_for_rgb(rgb: RGB): RGB {
  let lab = linear_srgb_to_oklab(rgb);

  let L = lab.L;
  let a = lab.a;
  let b = lab.b;

  let h = Math.atan2(b, a);
  let C = maximum_chroma_for_lh(L, h);

  return oklab_to_linear_srgb({ L: L, a: C * Math.cos(h), b: C * Math.sin(h) });
}

function print_color(c: RGB): void {
  console.log(linear_rgb_as_hex(c));
  let lab = linear_srgb_to_oklab(c);
  let C = Math.sqrt(lab.a * lab.a + lab.b * lab.b);
  let h = Math.atan2(lab.b, lab.a);
  console.log(`oklch(${lab.L}, ${C}, ${h})`);
}

function parse_hexcode(hex: string): RGB {
  let r = parseInt(hex.slice(1, 3), 16) / 255;
  let g = parseInt(hex.slice(3, 5), 16) / 255;
  let b = parseInt(hex.slice(5, 7), 16) / 255;

  return srgb_to_linear_rgb({ r: r, g: g, b: b });
}

module.exports = {
    linear_srgb_to_oklab,
    oklab_to_linear_srgb,
    srgb_to_linear_rgb,
    linear_rgb_to_srgb,
    rgb_as_hex,
    linear_rgb_as_hex,
    compute_max_saturation,
    find_cusp,
    find_gamut_intersection,
    clamp,
    sgn,
    gamut_clip_preserve_chroma,
    maximum_chroma_for_lh,
    maximize_chroma_for_rgb,
    print_color,
    parse_hexcode,
};

/*
>>> const X = require("./chroma");
>>> Z.print_color(Z.maximize_chroma_for_rgb(Z.parse_hexcode("#583e61")))
#72008d
oklch(0.40698242333583606, 0.19808518019193927, -0.7473789762405658)
*/

Gist

There's a super short example at the end that should illustrate how it works. In your specific case you probably want to call maximum_chroma_for_lh directly.

Answer 2

It's not an elegant solution but if you can convert from Oklch to an RGB space, you can increment the chroma until the conversion results in an RGB that contains a value outside of the 0-1 range.

For example, with my Unicolour library I can do this in C#:

const double increment = 0.001;

var chroma = 0.0;
while (true)
{
    // by default assumes sRGB display
    var unicolour = new Unicolour(ColourSpace.Oklch, 0.5, chroma, 285);
    if (!unicolour.IsInDisplayGamut)
    {
        chroma -= increment;
        break;
    }
    
    chroma += increment;
}

Console.WriteLine(chroma); // 0.292