Basic implementation of toGamut() (#199)

Author: danburzo

docs/api.md

@@ -376,7 +376,7 @@ By default, the color is converted to `lch` to perform the clamping, but any col
 import { clampChroma } from 'culori';
 
 clampChroma({ mode: 'oklch', l: 0.5, c: 0.16, h: 180 }, 'oklch');
-// => { mode: 'oklch', l: 0.5, c: 0.09, h: 180 }
+// ⇒ { mode: 'oklch', l: 0.5, c: 0.09, h: 180 }
 ```
 
 In general, chroma clamping is more accurate and computationally simpler when performed in the color's original space, where possible. Here's some sample code that uses the color's own `mode` for color spaces containing a Chroma dimension, and `lch` otherwise:
@@ -391,6 +391,31 @@ If the chroma-finding algorithm fails to find a displayable color (which can hap
 
 The function uses [the bisection method](https://en.wikipedia.org/wiki/Bisection_method) to speed up the search for the largest Chroma value. However, due to discontinuities in the CIELCh color space, the function is not guaranteed to return the optimal result. [See this discussion](https://github.com/d3/d3-color/issues/33) for details.
 
+<a id="toGamut" href="#toGamut">#</a> **toGamut**(_dest = 'rgb'_, _mode = 'oklch'_, _delta = differenceEuclidean('oklch')_, _jnd = 0.02_) → _function (color | string)_
+
+<span aria-label='Source:'>☞</span> [src/clamp.js]({{codebase}}/src/clamp.js)
+
+Obtain a color that's in the `dest` gamut, by first converting it to the `mode` color space and then finding the largest chroma that's in gamut, similar to `clampChroma()`.
+
+The color returned is in the `dest` color space.
+
+```js
+import { p3, toGamut } from 'culori';
+
+const color = 'lch(80% 150 60)';
+
+p3(color);
+// ⇒ { mode: "p3", r: 1.229…, g: 0.547…, b: -0.073… }
+
+const toP3 = toGamut('p3');
+toP3(color);
+// ⇒ { mode: "p3", r: 0.999…, g: 0.696…, b: 0.508… }
+```
+
+To address the shortcomings of `clampChroma`, which can sometimes produce colors more desaturated than necessary, the test used in the binary search is replaced with "is color is roughly in gamut", by comparing the candidate to the clipped version (obtained with `clampGamut`). The test passes if the colors are not to dissimilar, judged by the `delta` color difference function and an associated `jnd` just-noticeable difference value.
+
+The default arguments for this function correspond to [the gamut mapping algorithm](https://drafts.csswg.org/css-color/#css-gamut-mapping) defined in the CSS Color Module Level 4 spec, but the algorithm itself is slightly different.
+
 ## Interpolation
 
 In any color space, colors occupy positions given by their values for each channel. Interpolating colors means tracing a line through the coordinates of these colors, and figuring out what colors reside on the line at various positions.

package.json

@@ -58,6 +58,7 @@
 	"scripts": {
 		"prepare": "git config core.hooksPath .git-hooks",
 		"test": "tape 'test/**/*.test.js' | tap-spec",
+		"start": "npx esbuild --servedir=.",
 		"build": "node build.js",
 		"benchmark": "node benchmark/index.js",
 		"prepublishOnly": "npm run lint && npm run build && npm run test",

src/clamp.js

@@ -5,10 +5,16 @@ import { differenceEuclidean } from './difference.js';
 
 const rgb = converter('rgb');
 const fixup_rgb = c => {
-	c.r = Math.max(0, Math.min(c.r, 1));
-	c.g = Math.max(0, Math.min(c.g, 1));
-	c.b = Math.max(0, Math.min(c.b, 1));
-	return c;
+	const res = {
+		mode: c.mode,
+		r: Math.max(0, Math.min(c.r, 1)),
+		g: Math.max(0, Math.min(c.g, 1)),
+		b: Math.max(0, Math.min(c.b, 1))
+	};
+	if (c.alpha !== undefined) {
+		res.alpha = c.alpha;
+	}
+	return res;
 };
 
 const inrange_rgb = c => {
@@ -153,3 +159,89 @@ export function clampChroma(color, mode = 'lch') {
 		displayable(clamped) ? clamped : { ...clamped, c: _last_good_c }
 	);
 }
+
+/*
+	Obtain a color that's in the `dest` gamut,
+	by first converting it to the `mode` color space
+	and then finding the largest chroma that's in gamut,
+	similar to `clampChroma`. 
+
+	The color returned is in the `dest` color space.
+
+	To address the shortcomings of `clampChroma`, which can
+	sometimes produce colors more desaturated than necessary,
+	the test used in the binary search is replaced with
+	"is color is roughly in gamut", by comparing the candidate 
+	to the clipped version (obtained with `clampGamut`).
+	The test passes if the colors are not to dissimilar, 
+	judged by the `delta` color difference function 
+	and an associated `jnd` just-noticeable difference value.
+
+	The default arguments for this function correspond to the
+	gamut mapping algorithm defined in CSS Color Level 4:
+	https://drafts.csswg.org/css-color/#css-gamut-mapping
+ */
+export function toGamut(
+	dest = 'rgb',
+	mode = 'oklch',
+	delta = differenceEuclidean('oklch'),
+	jnd = 0.02
+) {
+	const destConv = converter(dest);
+
+	if (!getMode(dest).gamut) {
+		return color => destConv(color);
+	}
+
+	const inDestinationGamut = inGamut(dest);
+	const clipToGamut = clampGamut(dest);
+
+	const ucs = converter(mode);
+	const { ranges } = getMode(mode);
+
+	const gamutDef = getMode(dest);
+	const White = destConv('white');
+	const Black = destConv('black');
+
+	return color => {
+		color = prepare(color);
+		if (color === undefined) {
+			return undefined;
+		}
+		const candidate = { ...ucs(color) };
+		if (candidate.l >= ranges.l[1]) {
+			const res = { ...White };
+			if (color.alpha !== undefined) {
+				res.alpha = color.alpha;
+			}
+			return res;
+		}
+		if (candidate.l <= ranges.l[0]) {
+			const res = { ...Black };
+			if (color.alpha !== undefined) {
+				res.alpha = color.alpha;
+			}
+			return res;
+		}
+		if (inDestinationGamut(candidate)) {
+			return destConv(candidate);
+		}
+		let start = 0;
+		let end = candidate.c;
+		let ε = (ranges.c[1] - ranges.c[0]) / 4000; // 0.0001 for oklch()
+		let clipped = clipToGamut(candidate);
+		while (end - start > ε) {
+			candidate.c = (start + end) * 0.5;
+			clipped = clipToGamut(candidate);
+			if (
+				inDestinationGamut(candidate) ||
+				(jnd > 0 && delta(candidate, clipped) <= jnd)
+			) {
+				start = candidate.c;
+			} else {
+				end = candidate.c;
+			}
+		}
+		return destConv(inDestinationGamut(candidate) ? candidate : clipped);
+	};
+}

src/index-fn.js

@@ -100,7 +100,8 @@ export {
 	inGamut,
 	clampRgb,
 	clampChroma,
-	clampGamut
+	clampGamut,
+	toGamut
 } from './clamp.js';
 export { default as nearest } from './nearest.js';
 export { useMode, getMode, useParser, removeParser } from './modes.js';

src/index.js

@@ -101,7 +101,8 @@ export {
 	inGamut,
 	clampRgb,
 	clampChroma,
-	clampGamut
+	clampGamut,
+	toGamut
 } from './clamp.js';
 export { default as nearest } from './nearest.js';
 export { useMode, getMode, useParser, removeParser } from './modes.js';

test/api.test.js

@@ -210,6 +210,7 @@ const API_FULL = [
 	'serializeHex8',
 	'serializeHsl',
 	'serializeRgb',
+	'toGamut',
 	'unlerp',
 	'useMode',
 	'useParser',
@@ -449,6 +450,7 @@ const API_FN = [
 	'serializeHex8',
 	'serializeHsl',
 	'serializeRgb',
+	'toGamut',
 	'unlerp',
 	'useMode',
 	'useParser',

test/clamp.test.js

@@ -4,7 +4,9 @@ import {
 	displayable,
 	inGamut,
 	clampGamut,
-	formatCss
+	formatCss,
+	differenceEuclidean,
+	toGamut
 } from '../src/index.js';
 
 tape('RGB', function (test) {
@@ -134,3 +136,27 @@ tape('clampGamut()', t => {
 
 	t.end();
 });
+
+tape('toGamut()', t => {
+	t.equal(
+		formatCss(toGamut('rgb')({ mode: 'oklch', l: 1.5, c: 0.2, h: 180 })),
+		'color(srgb 1 1 1)',
+		'white'
+	);
+	t.equal(
+		formatCss(toGamut('rgb')({ mode: 'oklch', l: -1.5, c: 0.2, h: 180 })),
+		'color(srgb 0 0 0)',
+		'black'
+	);
+	t.equal(
+		formatCss(toGamut('rgb')('color(--lch-d65 100 0 180)')),
+		'color(srgb 0.9999999999999968 1.0000000000000016 0.9999999999999986)',
+		'chroma = 0'
+	);
+	t.equal(
+		formatCss(toGamut('p3')('lch(80% 150 60)')),
+		'color(display-p3 0.9999999999999994 0.6969234154991887 0.5084794582132421)',
+		'api docs example'
+	);
+	t.end();
+});

test/visual/gamut-map-newton.js

@@ -0,0 +1,89 @@
+const OKLAB_TO_LMS = [
+	[0.99999999845051981432, 0.39633779217376785678, 0.21580375806075880339],
+	[1.0000000088817607767, -0.1055613423236563494, -0.063854174771705903402],
+	[1.0000000546724109177, -0.089484182094965759684, -1.2914855378640917399]
+];
+
+const LMS_TO_XYZ_D65 = [
+	[1.2268798733741557, -0.5578149965554813, 0.28139105017721583],
+	[-0.04057576262431372, 1.1122868293970594, -0.07171106666151701],
+	[-0.07637294974672142, -0.4214933239627914, 1.5869240244272418]
+];
+
+const XYZ_D65_TO_LRGB = [
+	[3.2409699419045226, -1.537383177570094, -0.4986107602930034],
+	[-0.9692436362808796, 1.8759675015077202, 0.04155505740717559],
+	[0.05563007969699366, -0.20397695888897652, 1.0569715142428786]
+];
+
+const LMS_TO_LRGB = multiplyMatrices(XYZ_D65_TO_LRGB, LMS_TO_XYZ_D65);
+
+function multiplyMatrices(a, b) {
+	let res = [[], [], []];
+	for (let i = 0; i < 3; i++) {
+		for (let j = 0; j < 3; j++) {
+			res[i][j] =
+				a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
+		}
+	}
+	return res;
+}
+
+function multiplyMatrixWithVector(m, v) {
+	return [
+		m[0][0] * v[0] + m[0][1] * v[1] + m[0][2] * v[2],
+		m[1][0] * v[0] + m[1][1] * v[1] + m[1][2] * v[2],
+		m[2][0] * v[0] + m[2][1] * v[1] + m[2][2] * v[2]
+	];
+}
+
+function multiplyComponentWise(arr1, arr2) {
+	return arr1.map((it, i) => it * arr2[i]);
+}
+
+function multiplyVectors(arr1, arr2) {
+	return arr1.reduce((acc, curr, i) => {
+		return acc + curr * arr2[i];
+	}, 0);
+}
+
+export function GamutMapNewton(original) {
+	const zero_a_b = [0, original[1], original[2]];
+	let alpha = 1;
+	let res_oklab = original.slice();
+	let res_rgb = multiplyMatrixWithVector(
+		LMS_TO_LRGB,
+		multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => it ** 3)
+	);
+	for (let comp = 0; comp < 3; comp++) {
+		if (res_rgb[comp] >= 0 && res_rgb[comp] <= 1) continue;
+		let target = res_rgb[comp] > 1 ? 1 : 0;
+		for (let iter = 0; iter < 6; iter++) {
+			let residual =
+				multiplyVectors(
+					LMS_TO_LRGB[comp],
+					multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(
+						it => it ** 3
+					)
+				) - target;
+			if (Math.abs(residual) < 1e-15) break;
+			let gradient = multiplyVectors(
+				LMS_TO_LRGB[comp],
+				multiplyComponentWise(
+					multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(
+						it => 3 * it ** 2
+					),
+					multiplyMatrixWithVector(OKLAB_TO_LMS, zero_a_b)
+				)
+			);
+			alpha -= residual / gradient;
+			res_oklab[1] = alpha * original[1];
+			res_oklab[2] = alpha * original[2];
+		}
+		res_rgb = multiplyMatrixWithVector(
+			LMS_TO_LRGB,
+			multiplyMatrixWithVector(OKLAB_TO_LMS, res_oklab).map(it => it ** 3)
+		);
+	}
+	return [res_oklab, res_rgb];
+}