Maximize chroma in Oklch without changing L and h

chroma.ts

// 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)
*/