README.md

Fast OKLab Gamut Mapping in JavaScript

Using OKLCH and sRGB. This also works for Display P3, however you'll have to pass the correct coefficients to the cusp & intersection utilities, and instead of using linear sRGB, you would convert to linear Display P3.

If you wish to support both sRGB and Display P3 gamut mapping, it is possible with a single function and a matrix parameter, as they both use the same coordinate range (0..1) and gamma transfer function.

import { OKLCH_to_linear_sRGB } from './convert.js';
import { rgb_in_gamut, linear_to_gamma } from './util.js';
import { find_cusp, find_gamut_intersection } from './oklab-cusp.js';

let lch = [ 0.5, 0.2, 270 ];
let lrgb = OKLCH_to_linear_sRGB(lch);

// if not exactly in gamut
if (!rgb_in_gamut(lrgb, 0)) {
  // do gamut mapping
  const LC = mapToGamutLcusp(lch, cusp);
  // get new LCH
  lch = [LC[0], LC[1], H];
  // convert again to lRGB
  lrgb = OKLCH_to_linear_sRGB(lch);
}

// note: after gamut mapping, it may be on the edge, technically outside of sRGB gamut
// but it is so close that we just clip it, as it is far below
// the "JND" in terms of DeltaEOK
const [ R, G, B ] = lrgb.map(n => clamp(linear_to_gamma(n), 0, 1));

// gamma sRGB coords in 0 ... 1
console.log(R,G,B);

function mapToGamutLcusp(
  LCH,
  cusp
) {
  const [L, C, H] = LCH;
  const hueAngle = degToRad(H);
  const aNorm = Math.cos(hueAngle);
  const bNorm = Math.sin(hueAngle);
  cusp = cusp || find_cusp(aNorm, bNorm);
    
  // Cusp luminance as target
  const LTarget = cusp[0];
  
  // Other strategies
  // const LTarget = clamp(LCH[0], 0, 1);
  // const LTarget = 0.5;

  const t = find_gamut_intersection(
    aNorm,
    bNorm,
    L,
    C,
    LTarget,
    cusp
  );
  const L_clipped = LTarget * (1 - t) + t * L;
  const C_clipped = t * C;
  return [L_clipped, C_clipped, H];
}

function clamp(value, min, max) {
  return Math.min(max, Math.max(min, value));
}

Credits

The matrices and conversion code has been mostly pulled from Colorjs and Coloraide. It exists mostly for reference, it is probably wiser to use the Colorjs library rather than relying on my conversion code here with no tests/coverage.

Code and reference for sRGB gamut clipping is from Björn Ottosson's blog on the subject here and here.

License

MIT License.

convert.js

import { mat3, vec3 } from "gl-matrix";
import {
  clamp,
  gamma_to_linear,
  linear_to_gamma,
  multiplyMatrices,
} from "./util.js";

// OKLab and OKLCH
// https://bottosson.github.io/posts/oklab/

// XYZ <-> LMS matrices recalculated for consistent reference white
// see https://github.com/w3c/csswg-drafts/issues/6642#issuecomment-943521484
// recalculated for 64bit precision
// see https://github.com/color-js/color.js/pull/357

// Given XYZ relative to D65, convert to OKLab
export const XYZ_to_LMS_M = [
  [0.819022437996703, 0.3619062600528904, -0.1288737815209879],
  [0.0329836539323885, 0.9292868615863434, 0.0361446663506424],
  [0.0481771893596242, 0.2642395317527308, 0.6335478284694309],
];

export const LMS_to_OKLab_M = [
  [0.210454268309314, 0.7936177747023054, -0.0040720430116193],
  [1.9779985324311684, -2.4285922420485799, 0.450593709617411],
  [0.0259040424655478, 0.7827717124575296, -0.8086757549230774],
];

// Given OKLab, convert to XYZ relative to D65
export const LMS_to_XYZ_M = [
  [1.2268798758459243, -0.5578149944602171, 0.2813910456659647],
  [-0.0405757452148008, 1.112286803280317, -0.0717110580655164],
  [-0.0763729366746601, -0.4214933324022432, 1.5869240198367816],
];

export const OKLab_to_LMS_M = [
  [1.0, 0.3963377773761749, 0.2158037573099136],
  [1.0, -0.1055613458156586, -0.0638541728258133],
  [1.0, -0.0894841775298119, -1.2914855480194092],
];

export const linear_sRGB_to_LMS_M = [
  [0.4122214694707629, 0.5363325372617349, 0.051445993267502196],
  [0.2119034958178251, 0.6806995506452345, 0.10739695353694051],
  [0.08830245919005637, 0.2817188391361215, 0.6299787016738223],
];

export const LMS_to_linear_sRGB_M = [
  [4.076741636075959, -3.307711539258062, 0.2309699031821041],
  [-1.2684379732850313, 2.6097573492876878, -0.3413193760026569],
  [-0.004196076138675526, -0.703418617935936, 1.7076146940746113],
];

export const LMS_to_linear_P3_M = [
  [3.127768971361874, -2.2571357625916395, 0.12936679122976516],
  [-1.0910090184377979, 2.413331710306922, -0.32232269186912466],
  [-0.02601080193857028, -0.508041331704167, 1.5340521336427373],
];

export const linear_P3_to_LMS_M = [
  [0.4813798527499543, 0.4621183710113182, 0.05650177623872754],
  [0.2288319418112447, 0.6532168193835677, 0.11795123880518772],
  [0.08394575232299314, 0.22416527097756647, 0.6918889766994405],
];

// https://github.com/w3c/csswg-drafts/issues/5922
export const linear_sRGB_to_XYZ_M = [
  [0.41239079926595934, 0.357584339383878, 0.1804807884018343],
  [0.21263900587151027, 0.715168678767756, 0.07219231536073371],
  [0.01933081871559182, 0.11919477979462598, 0.9505321522496607],
];

export const XYZ_to_linear_sRGB_M = [
  [3.2409699419045226, -1.537383177570094, -0.4986107602930034],
  [-0.9692436362808796, 1.8759675015077202, 0.04155505740717559],
  [0.05563007969699366, -0.20397695888897652, 1.0569715142428786],
];

// Display-P3 to XYZ and back
export const linear_P3_to_XYZ_M = [
  [608311 / 1250200, 189793 / 714400, 198249 / 1000160],
  [35783 / 156275, 247089 / 357200, 198249 / 2500400],
  [0 / 1, 32229 / 714400, 5220557 / 5000800],
];

export const XYZ_to_linear_P3_M = [
  [446124 / 178915, -333277 / 357830, -72051 / 178915],
  [-14852 / 17905, 63121 / 35810, 423 / 17905],
  [11844 / 330415, -50337 / 660830, 316169 / 330415],
];

const floatToByte = (n) => clamp(Math.round(255 * n), 0, 255);

export function XYZ_to_OKLab(XYZ) {
  const LMS = multiplyMatrices(XYZ_to_LMS_M, XYZ);
  // JavaScript Math.cbrt returns a sign-matched cube root
  // beware if porting to other languages
  // especially if tempted to use a general power function
  return multiplyMatrices(
    LMS_to_OKLab_M,
    LMS.map((c) => Math.cbrt(c))
  );
  // L in range [0,1]. For use in CSS, multiply by 100 and add a percent
}

export function OKLab_to_XYZ(OKLab) {
  const LMSnl = multiplyMatrices(OKLab_to_LMS_M, OKLab);
  return multiplyMatrices(
    LMS_to_XYZ_M,
    LMSnl.map((c) => c ** 3)
  );
}

// L = 0...1
// a, b = -0.4...0.4
export function OKLab_to_OKLCH(OKLab) {
  const hue = (Math.atan2(OKLab[2], OKLab[1]) * 180) / Math.PI;
  return [
    OKLab[0], // L is still L
    Math.sqrt(OKLab[1] ** 2 + OKLab[2] ** 2), // Chroma
    hue >= 0 ? hue : hue + 360, // Hue, in degrees [0 to 360)
  ];
}

// L = 0...1
// C = -0.4..0.4
// H = 0...360
export function OKLCH_to_OKLab(OKLCH) {
  return [
    OKLCH[0], // L is still L
    OKLCH[1] * Math.cos((OKLCH[2] * Math.PI) / 180), // a
    OKLCH[1] * Math.sin((OKLCH[2] * Math.PI) / 180), // b
  ];
}

export function linear_sRGB_to_XYZ(sRGB) {
  // convert an array of linear-light sRGB values to CIE XYZ
  // using sRGB's own white, D65 (no chromatic adaptation)
  return multiplyMatrices(linear_sRGB_to_XYZ_M, sRGB);
}

export function XYZ_to_linear_sRGB(XYZ) {
  // convert XYZ to linear-light sRGB
  return multiplyMatrices(XYZ_to_linear_sRGB_M, XYZ);
}

export const linear_sRGB_to_sRGB = (sRGB) =>
  sRGB.map((n) => floatToByte(linear_to_gamma(n)));

export const sRGB_to_linear_sRGB = (sRGB) =>
  sRGB.map((n) => gamma_to_linear(n / 0xff));

export function linear_P3_to_XYZ(p3) {
  // convert an array of linear-light display-p3 values to CIE XYZ
  // using  D65 (no chromatic adaptation)
  // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
  return multiplyMatrices(linear_P3_to_XYZ_M, p3);
}

export function XYZ_to_linear_P3(XYZ) {
  // convert XYZ to linear-light P3
  return multiplyMatrices(XYZ_to_linear_P3_M, XYZ);
}

export const linear_P3_to_P3 = (p3) => p3.map((n) => linear_to_gamma(n));
export const P3_to_linear_P3 = (p3) => p3.map((n) => gamma_to_linear(n));

// const spaces = [
//   'srgb',
//   'srgb-linear',
//   'xyz',
//   'xyy',
//   'display-p3',
//   'oklab',
//   'oklch'
// ];

// function convert (fromCoords, fromId, toId) {
// }

////// Utilities
const epsilon = 0.000075;

export const linear_sRGB_in_gamut = (lrgb) =>
  lrgb.every((n) => {
    n = linear_to_gamma(n);
    return n >= -epsilon && n <= 1 + epsilon;
  });

export function linear_sRGB_to_OKLab(sRGB) {
  const R = sRGB[0];
  const G = sRGB[1];
  const B = sRGB[2];
  const [l, m, s] = linear_sRGB_to_LMS_M.map((vec) => {
    return vec[0] * R + vec[1] * G + vec[2] * B;
  });
  const l_ = Math.cbrt(l);
  const m_ = Math.cbrt(m);
  const s_ = Math.cbrt(s);
  return LMS_to_OKLab_M.map((vec) => {
    return vec[0] * l_ + vec[1] * m_ + vec[2] * s_;
  });
}

export function OKLab_to_linear_sRGB(OKLab) {
  const L = OKLab[0];
  const A = OKLab[1];
  const B = OKLab[2];
  const [l_, m_, s_] = OKLab_to_LMS_M.map((vec) => {
    return vec[0] * L + vec[1] * A + vec[2] * B;
  });

  const l = l_ * l_ * l_;
  const m = m_ * m_ * m_;
  const s = s_ * s_ * s_;
  return LMS_to_linear_sRGB_M.map((vec) => {
    return vec[0] * l + vec[1] * m + vec[2] * s;
  });
}

export function OKLab_to_linear_P3(OKLab) {
  const L = OKLab[0];
  const A = OKLab[1];
  const B = OKLab[2];
  const [l_, m_, s_] = OKLab_to_LMS_M.map((vec) => {
    return vec[0] * L + vec[1] * A + vec[2] * B;
  });

  const l = l_ * l_ * l_;
  const m = m_ * m_ * m_;
  const s = s_ * s_ * s_;
  return LMS_to_linear_P3_M.map((vec) => {
    return vec[0] * l + vec[1] * m + vec[2] * s;
  });
}

export function linear_P3_to_OKLab(P3) {
  const R = P3[0];
  const G = P3[1];
  const B = P3[2];
  const [l, m, s] = linear_P3_to_LMS_M.map((vec) => {
    return vec[0] * R + vec[1] * G + vec[2] * B;
  });
  const l_ = Math.cbrt(l);
  const m_ = Math.cbrt(m);
  const s_ = Math.cbrt(s);
  return LMS_to_OKLab_M.map((vec) => {
    return vec[0] * l_ + vec[1] * m_ + vec[2] * s_;
  });
}

export const OKLCH_to_linear_sRGB = (OKLCH) =>
  OKLab_to_linear_sRGB(OKLCH_to_OKLab(OKLCH));

export const OKLab_to_sRGB = (OKLab) =>
  linear_sRGB_to_sRGB(OKLab_to_linear_sRGB(OKLab));

export const OKLCH_to_sRGB = (OKLCH) =>
  linear_sRGB_to_sRGB(OKLCH_to_linear_sRGB(OKLCH));

// export const linear_sRGB_to_OKLab = (sRGB) =>
//   XYZ_to_OKLab(linear_sRGB_to_XYZ(sRGB));

export const linear_sRGB_to_OKLCH = (sRGB) =>
  OKLab_to_OKLCH(linear_sRGB_to_OKLab(sRGB));

export const sRGB_to_OKLab = (sRGB) =>
  linear_sRGB_to_OKLab(sRGB_to_linear_sRGB(sRGB));

export const sRGB_to_OKLCH = (sRGB) =>
  linear_sRGB_to_OKLCH(sRGB_to_linear_sRGB(sRGB));

// p3 - oklab

// direct versions now exist
// export const OKLab_to_linear_P3 = (OKLab) =>
//   XYZ_to_linear_P3(OKLab_to_XYZ(OKLab));
// export const linear_P3_to_OKLab = (p3) => XYZ_to_OKLab(linear_P3_to_XYZ(p3));

export const OKLab_to_P3 = (OKLab) =>
  linear_P3_to_P3(OKLab_to_linear_P3(OKLab));

export const OKLCH_to_linear_P3 = (OKLCH) =>
  XYZ_to_linear_P3(OKLab_to_XYZ(OKLCH_to_OKLab(OKLCH)));

export const OKLCH_to_P3 = (OKLCH) =>
  linear_P3_to_P3(OKLCH_to_linear_P3(OKLCH));

export const P3_to_OKLab = (p3) => linear_P3_to_OKLab(P3_to_linear_P3(p3));

export const XYZ_to_sRGB = (XYZ) =>
  linear_sRGB_to_sRGB(XYZ_to_linear_sRGB(XYZ));

oklab-cusp.js

//// OKHSL and OKHSV
import {
  OKLab_to_linear_sRGB,
  OKLab_to_LMS_M,
  LMS_to_linear_sRGB_M,
} from "./convert.js";

import { vec2 } from "gl-matrix";

export const sRGB_coefficients = [
  // Red
  [
    // Limit
    [-1.8817030993265886, -0.8093650129914305],
    // `Kn` coefficients
    [1.19086277, 1.76576728, 0.59662641, 0.75515197, 0.56771245],
  ],
  // Green
  [
    // Limit
    [1.8144407988011002, -1.1944526678052347],
    // `Kn` coefficients
    [0.73956515, -0.45954404, 0.08285427, 0.12541073, -0.14503204],
  ],
  // Blue
  [
    // Limit
    [0.13110757611180976, 1.8133397092666066],
    // `Kn` coefficients
    [1.35733652, -0.00915799, -1.1513021, -0.50559606, 0.00692167],
  ],
];

export const P3_coefficients = [
  // Red
  [
    // Limit R_DIR
    [-1.7723439275129815, -0.8207587433674075],
    // `Kn` coefficients R
    [
      1.1941401817279191, 1.7629811997119313, 0.5958599382482035,
      0.7575999740543717, 0.5681684967815825,
    ],
  ],
  // Green
  [
    // Limit G_DIR
    [1.8031987175305486, -1.1932813966558908],
    // `Kn` coefficients G
    [
      0.739566819222909, -0.45954279991708363, 0.08285308769183684,
      0.12541164951759598, -0.14503290744198968,
    ],
  ],
  // Blue
  [
    // Limit B_DIR
    [0.08970487824467556, 1.903277465741611],
    // `Kn` coefficients B
    [
      1.365094411769746, -0.013962295570417609, -1.145230508988296,
      -0.5025987876722386, 0.003174713115376433,
    ],
  ],
];

const floatMax = Number.MAX_VALUE;
const K1 = 0.206;
const K2 = 0.03;
const K3 = (1.0 + K1) / (1.0 + K2);

export function compute_max_saturation(
  a,
  b,
  lmsToRgb = LMS_to_linear_sRGB_M,
  okCoeff = sRGB_coefficients
) {
  // https://github.com/color-js/color.js/blob/main/src/spaces/okhsl.js
  // Finds the maximum saturation possible for a given hue that fits in RGB.
  //
  // Saturation here is defined as `S = C/L`.
  // `a` and `b` must be normalized so `a^2 + b^2 == 1`.

  // Max saturation will be when one of r, g or b goes below zero.

  // Select different coefficients depending on which component goes below zero first.

  let k0, k1, k2, k3, k4, wl, wm, ws;

  if (vec2.dot(okCoeff[0][0], [a, b]) > 1) {
    // Red component
    [k0, k1, k2, k3, k4] = okCoeff[0][1];
    [wl, wm, ws] = lmsToRgb[0];
  } else if (vec2.dot(okCoeff[1][0], [a, b]) > 1) {
    // Green component
    [k0, k1, k2, k3, k4] = okCoeff[1][1];
    [wl, wm, ws] = lmsToRgb[1];
  } else {
    // Blue component
    [k0, k1, k2, k3, k4] = okCoeff[2][1];
    [wl, wm, ws] = lmsToRgb[2];
  }

  // Approximate max saturation using a polynomial:
  let sat = k0 + k1 * a + k2 * b + k3 * a ** 2 + k4 * a * b;

  // Do one step Halley's method to get closer.
  // This gives an error less than 10e6, except for some blue hues where the `dS/dh` is close to infinite.
  // This should be sufficient for most applications, otherwise do two/three steps.

  let kl = vec2.dot(OKLab_to_LMS_M[0].slice(1), [a, b]);
  let km = vec2.dot(OKLab_to_LMS_M[1].slice(1), [a, b]);
  let ks = vec2.dot(OKLab_to_LMS_M[2].slice(1), [a, b]);

  let l_ = 1.0 + sat * kl;
  let m_ = 1.0 + sat * km;
  let s_ = 1.0 + sat * ks;

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

  let lds = 3.0 * kl * l_ ** 2;
  let mds = 3.0 * km * m_ ** 2;
  let sds = 3.0 * ks * s_ ** 2;

  let lds2 = 6.0 * kl ** 2 * l_;
  let mds2 = 6.0 * km ** 2 * m_;
  let sds2 = 6.0 * ks ** 2 * s_;

  let f = wl * l + wm * m + ws * s;
  let f1 = wl * lds + wm * mds + ws * sds;
  let f2 = wl * lds2 + wm * mds2 + ws * sds2;

  sat = sat - (f * f1) / (f1 ** 2 - 0.5 * f * f2);

  return sat;
}

export function find_cusp(
  a,
  b,
  lmsToRgb,
  okCoeff,
  OKLab_to_linear_RGB = OKLab_to_linear_sRGB
) {
  // First, find the maximum saturation (saturation S = C/L)
  var S_cusp = compute_max_saturation(a, b, lmsToRgb, okCoeff);
  // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1:
  var rgb_at_max = OKLab_to_linear_RGB([1, S_cusp * a, S_cusp * b]);
  var L_cusp = Math.cbrt(
    1 / Math.max(Math.max(rgb_at_max[0], rgb_at_max[1]), rgb_at_max[2])
  );
  var C_cusp = L_cusp * S_cusp;
  return [L_cusp, C_cusp];
}

export function find_gamut_intersection(
  a,
  b,
  l1,
  c1,
  l0,
  cusp,
  lmsToRgb,
  okCoeff
) {
  // Finds intersection of the line.
  //
  // Defined by the following:
  //
  // ```
  // L = L0 * (1 - t) + t * L1
  // C = t * C1
  // ```
  //
  // `a` and `b` must be normalized so `a^2 + b^2 == 1`.

  let t;

  if (!lmsToRgb) throw new Error("Must pass lmsToRgb");

  if (!cusp) {
    cusp = find_cusp(a, b, lmsToRgb, okCoeff);
  }

  // Find the intersection for upper and lower half separately
  if ((l1 - l0) * cusp[1] - (cusp[0] - l0) * c1 <= 0.0) {
    // Lower half
    t = (cusp[1] * l0) / (c1 * cusp[0] + cusp[1] * (l0 - l1));
  } else {
    // Upper half

    // First intersect with triangle
    t = (cusp[1] * (l0 - 1.0)) / (c1 * (cusp[0] - 1.0) + cusp[1] * (l0 - l1));

    // Then one step Halley's method
    let dl = l1 - l0;
    let dc = c1;

    let kl = vec2.dot(OKLab_to_LMS_M[0].slice(1), [a, b]);
    let km = vec2.dot(OKLab_to_LMS_M[1].slice(1), [a, b]);
    let ks = vec2.dot(OKLab_to_LMS_M[2].slice(1), [a, b]);

    let ldt_ = dl + dc * kl;
    let mdt_ = dl + dc * km;
    let sdt_ = dl + dc * ks;

    // If higher accuracy is required, 2 or 3 iterations of the following block can be used:
    let L = l0 * (1.0 - t) + t * l1;
    let C = t * c1;

    let l_ = L + C * kl;
    let m_ = L + C * km;
    let s_ = L + C * ks;

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

    let ldt = 3 * ldt_ * l_ ** 2;
    let mdt = 3 * mdt_ * m_ ** 2;
    let sdt = 3 * sdt_ * s_ ** 2;

    let ldt2 = 6 * ldt_ ** 2 * l_;
    let mdt2 = 6 * mdt_ ** 2 * m_;
    let sdt2 = 6 * sdt_ ** 2 * s_;

    let r_ = vec2.dot(lmsToRgb[0], [l, m, s]) - 1;
    let r1 = vec2.dot(lmsToRgb[0], [ldt, mdt, sdt]);
    let r2 = vec2.dot(lmsToRgb[0], [ldt2, mdt2, sdt2]);

    let ur = r1 / (r1 * r1 - 0.5 * r_ * r2);
    let tr = -r_ * ur;

    let g_ = vec2.dot(lmsToRgb[1], [l, m, s]) - 1;
    let g1 = vec2.dot(lmsToRgb[1], [ldt, mdt, sdt]);
    let g2 = vec2.dot(lmsToRgb[1], [ldt2, mdt2, sdt2]);

    let ug = g1 / (g1 * g1 - 0.5 * g_ * g2);
    let tg = -g_ * ug;

    let b_ = vec2.dot(lmsToRgb[2], [l, m, s]) - 1;
    let b1 = vec2.dot(lmsToRgb[2], [ldt, mdt, sdt]);
    let b2 = vec2.dot(lmsToRgb[2], [ldt2, mdt2, sdt2]);

    let ub = b1 / (b1 * b1 - 0.5 * b_ * b2);
    let tb = -b_ * ub;

    tr = ur >= 0.0 ? tr : floatMax;
    tg = ug >= 0.0 ? tg : floatMax;
    tb = ub >= 0.0 ? tb : floatMax;

    t += Math.min(tr, Math.min(tg, tb));
  }

  return t;
}

export function get_ST_max(a_, b_, cusp) {
  if (cusp === void 0) {
    cusp = null;
  }
  if (!cusp) {
    cusp = find_cusp(a_, b_);
  }
  var L = cusp[0];
  var C = cusp[1];
  return [C / L, C / (1 - L)];
}

util.js

const roundByte = (n) => clamp(Math.round(n), 0, 255);

/**
 * A is m x n. B is n x p. product is m x p.
 * @param {number[] | number[][]} A Matrix m x n or a vector
 * @param {number[] | number[][]} B Matrix n x p or a vector
 * @returns {number[]} Matrix m x p
 */
export function multiplyMatrices(A, B) {
  let m = A.length;

  if (!Array.isArray(A[0])) {
    // A is vector, convert to [[a, b, c, ...]]
    A = [A];
  }

  if (!Array.isArray(B[0])) {
    // B is vector, convert to [[a], [b], [c], ...]]
    B = B.map((x) => [x]);
  }

  let p = B[0].length;
  let B_cols = B[0].map((_, i) => B.map((x) => x[i])); // transpose B
  let product = A.map((row) =>
    B_cols.map((col) => {
      let ret = 0;

      if (!Array.isArray(row)) {
        for (let c of col) {
          ret += row * c;
        }

        return ret;
      }

      for (let i = 0; i < row.length; i++) {
        ret += row[i] * (col[i] || 0);
      }

      return ret;
    })
  );

  if (m === 1) {
    product = product[0]; // Avoid [[a, b, c, ...]]
  }

  if (p === 1) {
    return product.map((x) => x[0]); // Avoid [[a], [b], [c], ...]]
  }

  return product;
}

export function gamma_to_linear(val) {
  // convert a single channel value
  // where in-gamut values are in the range [0 - 1]
  // to linear light (un-companded) form.
  // https://en.wikipedia.org/wiki/SRGB
  // Extended transfer function:
  // for negative values,  linear portion is extended on reflection of axis,
  // then reflected power function is used.
  let sign = val < 0 ? -1 : 1;
  let abs = Math.abs(val);
  if (abs <= 0.04045) {
    return val / 12.92;
  }
  return sign * Math.pow((abs + 0.055) / 1.055, 2.4);
}

export function linear_to_gamma(val) {
  // convert a single channel linear-light value in range 0-1
  // to gamma corrected form
  // https://en.wikipedia.org/wiki/SRGB
  // Extended transfer function:
  // For negative values, linear portion extends on reflection
  // of axis, then uses reflected pow below that
  let sign = val < 0 ? -1 : 1;
  let abs = Math.abs(val);
  if (abs > 0.0031308) {
    return sign * (1.055 * Math.pow(abs, 1 / 2.4) - 0.055);
  }
  return 12.92 * val;
}

export const clamp = (value, min, max) => Math.max(Math.min(value, max), min);

export const normalizeHue = (hue) => ((hue = hue % 360) < 0 ? hue + 360 : hue);

export const hex_to_rgb = (str) => {
  let hex = str.replace(/#/, "");
  if (hex.length === 3) {
    // expand shorthand
    hex = hex[0] + hex[0] + hex[1] + hex[1] + hex[2] + hex[2];
  } else if (hex.length > 6) {
    // discard alpha
    hex = hex.slice(0, 6);
  }
  const rgb = parseInt(hex, 16);
  let r = (rgb >> 16) & 0xff;
  let g = (rgb >> 8) & 0xff;
  let b = rgb & 0xff;
  return [r, g, b];
};

export const rgb_to_hex = (rgb) =>
  `#${rgb.map((n) => roundByte(n).toString(16).padStart(2, "0")).join("")}`;

// const GAMUT_EPSILON = 1e-6;
// const GAMUT_MIN = -GAMUT_EPSILON;
// const GAMUT_MAX = 1 + GAMUT_EPSILON;

const epsilon = 0.000075;

export const rgb_in_gamut = (lrgb, ep = epsilon) =>
  lrgb.every((n) => n >= -ep && n <= 1 + ep);

// in degrees
export const angle_delta = (angle1, angle2) => {
  const diff = ((angle2 - angle1 + 180) % 360) - 180;
  return diff < -180 ? diff + 360 : diff;
};

export const xyY_to_XYZ = (arg) => {
  var X, Y, Z, x, y;
  x = arg[0];
  y = arg[1];
  Y = arg[2];
  if (y === 0) {
    return [0, 0, 0];
  }
  X = (x * Y) / y;
  Z = ((1 - x - y) * Y) / y;
  return [X, Y, Z];
};

export const XYZ_to_xyY = (arg) => {
  var sum, X, Y, Z;
  X = arg[0];
  Y = arg[1];
  Z = arg[2];
  sum = X + Y + Z;
  if (sum === 0) {
    return [0, 0, Y];
  }
  return [X / sum, Y / sum, Y];
};